Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/master' into fluxDiagram_Chosen_…
Browse files Browse the repository at this point in the history 
…Species

Also, added declaration of global variables in methods for the global
variables like maximumNodeCount to avoid interpreting them as local
variables.
  • Loading branch information
Nick Vandewiele committed Jun 13, 2012
2 parents 70e0335 + 4fe5933 commit 7e65555
Show file tree
Hide file tree
Showing 2 changed files with 49 additions and 29 deletions.
76 changes: 47 additions & 29 deletions generateFluxDiagram.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,14 +47,22 @@

################################################################################

def generateFluxDiagram(reactionModel, times, concentrations, reactionRates, outputDirectory, centralSpecies = None, speciesDirectory=None):
def generateFluxDiagram(reactionModel, times, concentrations, reactionRates, outputDirectory, centralSpecies = None, speciesDirectory=None, settings = None):
"""
For a given `reactionModel` and simulation results stored as arrays of
`times`, species `concentrations`, and `reactionRates`, generate a series
of flux diagrams as frames of an animation, then stitch them together into
a movie. The individual frames and the final movie are saved on disk at
`outputDirectory.`
"""
global maximumNodeCount, maximumEdgeCount, timeStep, concentrationTolerance, speciesRateTolerance
# Allow user defined settings for flux diagram generation if given
if settings:
maximumNodeCount = settings['maximumNodeCount']
maximumEdgeCount = settings['maximumEdgeCount']
timeStep = settings['timeStep']
concentrationTolerance = settings['concentrationTolerance']
speciesRateTolerance = settings['speciesRateTolerance']

# Get the species and reactions corresponding to the provided concentrations and reaction rates
speciesList = reactionModel.core.species[:]
Expand Down Expand Up @@ -153,11 +161,12 @@ def generateFluxDiagram(reactionModel, times, concentrations, reactionRates, out
node.set_label("")
# Add an edge for each species-species rate
for reactantIndex, productIndex in edges:
reactant = speciesList[reactantIndex]
product = speciesList[productIndex]
edge = pydot.Edge(reactant.label, product.label)
edge.set_penwidth(maximumEdgePenWidth)
graph.add_edge(edge)
if reactantIndex in nodes and productIndex in nodes:
reactant = speciesList[reactantIndex]
product = speciesList[productIndex]
edge = pydot.Edge(reactant.label, product.label)
edge.set_penwidth(maximumEdgePenWidth)
graph.add_edge(edge)

# Generate the coordinates for all of the nodes using the specified program
graph = pydot.graph_from_dot_data(graph.create_dot(prog=program))
Expand All @@ -181,23 +190,24 @@ def generateFluxDiagram(reactionModel, times, concentrations, reactionRates, out
slope = -maximumEdgePenWidth / math.log10(speciesRateTolerance)
for index in range(len(edges)):
reactantIndex, productIndex = edges[index]
reactant = speciesList[reactantIndex]
product = speciesList[productIndex]
edge = graph.get_edge('"{0}"'.format(reactant.label), '"{0}"'.format(product.label))[0]
# Determine direction of arrow based on sign of rate
speciesRate = speciesRates[t,reactantIndex,productIndex] / maxSpeciesRate
if speciesRate < 0:
edge.set_dir("back")
speciesRate = -speciesRate
else:
edge.set_dir("forward")
# Set the edge pen width
if speciesRate < speciesRateTolerance:
penwidth = 0.0
edge.set_dir("none")
else:
penwidth = slope * math.log10(speciesRate) + maximumEdgePenWidth
edge.set_penwidth(penwidth)
if reactantIndex in nodes and productIndex in nodes:
reactant = speciesList[reactantIndex]
product = speciesList[productIndex]
edge = graph.get_edge('"{0}"'.format(reactant.label), '"{0}"'.format(product.label))[0]
# Determine direction of arrow based on sign of rate
speciesRate = speciesRates[t,reactantIndex,productIndex] / maxSpeciesRate
if speciesRate < 0:
edge.set_dir("back")
speciesRate = -speciesRate
else:
edge.set_dir("forward")
# Set the edge pen width
if speciesRate < speciesRateTolerance:
penwidth = 0.0
edge.set_dir("none")
else:
penwidth = slope * math.log10(speciesRate) + maximumEdgePenWidth
edge.set_penwidth(penwidth)
# Save the graph at this time to a dot file and a PNG image
if times[t] == 0:
label = 't = 0 s'
Expand Down Expand Up @@ -234,11 +244,19 @@ def generateFluxDiagram(reactionModel, times, concentrations, reactionRates, out

################################################################################

def simulate(reactionModel, reactionSystem):
def simulate(reactionModel, reactionSystem, settings = None):
"""
Generate and return a set of core and edge species and reaction fluxes
by simulating the given `reactionSystem` using the given `reactionModel`.
"""
global maximumNodeCount, maximumEdgeCount, timeStep, concentrationTolerance, speciesRateTolerance
# Allow user defined settings for flux diagram generation if given
if settings:
maximumNodeCount = settings['maximumNodeCount']
maximumEdgeCount = settings['maximumEdgeCount']
timeStep = settings['timeStep']
concentrationTolerance = settings['concentrationTolerance']
speciesRateTolerance = settings['speciesRateTolerance']

coreSpecies = reactionModel.core.species
coreReactions = reactionModel.core.reactions
Expand Down Expand Up @@ -494,7 +512,7 @@ def loadRMGPyJob(inputFile, chemkinFile=None, speciesDict=None):

################################################################################

def createFluxDiagram(savePath, inputFile, chemkinFile, speciesDict, java = False, chemkinOutput = '', centralSpecies = None):
def createFluxDiagram(savePath, inputFile, chemkinFile, speciesDict, java = False, settings = None, chemkinOutput = '', centralSpecies = None):
"""
Generates the flux diagram based on a condition 'inputFile', chemkin.inp chemkinFile,
a speciesDict txt file, plus an optional chemkinOutput file.
Expand All @@ -517,7 +535,7 @@ def createFluxDiagram(savePath, inputFile, chemkinFile, speciesDict, java = Fals
time, coreSpeciesConcentrations, coreReactionRates, edgeReactionRates = loadChemkinOutput(chemkinOutput, rmg.reactionModel)

print 'Generating flux diagram for chemkin output...'
generateFluxDiagram(rmg.reactionModel, time, coreSpeciesConcentrations, coreReactionRates, os.path.join(savePath, '1'), centralSpecies, speciesPath)
generateFluxDiagram(rmg.reactionModel, time, coreSpeciesConcentrations, coreReactionRates, os.path.join(savePath, '1'), centralSpecies, speciesPath, settings)

else:
# Generate a flux diagram video for each reaction system
Expand All @@ -536,11 +554,11 @@ def createFluxDiagram(savePath, inputFile, chemkinFile, speciesDict, java = Fals
reactionSystem.termination.append(TerminationTime((1e10,'s')))

print 'Conducting simulation of reaction system {0:d}...'.format(index+1)
time, coreSpeciesConcentrations, coreReactionRates, edgeReactionRates = simulate(rmg.reactionModel, reactionSystem)
time, coreSpeciesConcentrations, coreReactionRates, edgeReactionRates = simulate(rmg.reactionModel, reactionSystem, settings)

print 'Generating flux diagram for reaction system {0:d}...'.format(index+1)
generateFluxDiagram(rmg.reactionModel, time, coreSpeciesConcentrations, coreReactionRates, os.path.join(savePath, '{0:d}'.format(index+1)), centralSpecies, speciesPath)

generateFluxDiagram(rmg.reactionModel, time, coreSpeciesConcentrations, coreReactionRates, os.path.join(savePath, '{0:d}'.format(index+1)),
centralSpecies, speciesPath, settings)
################################################################################

if __name__ == '__main__':
Expand Down
2 changes: 2 additions & 0 deletions rmgpy/cantherm/gaussian.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,6 +146,8 @@ def loadGeometry(self):
mass[i] = 12.0
elif number[i] == 8:
mass[i] = 15.99491461956
elif number[i] == 16:
mass[i] = 31.97207100
else:
print 'Atomic number %i not yet supported in loadGeometry().' % number[i]

Expand Down

0 comments on commit 7e65555

Please sign in to comment.