Skip to content

Commit

Permalink
Merge pull request #364 from ReactionMechanismGenerator/rmg2to3scripts
Browse files Browse the repository at this point in the history 
Update RMG-database scripts and notebooks to Python 3
  • Loading branch information
@mliu49
mliu49 authored Dec 4, 2019
2 parents 3d0ea81 + 7f5fa5c commit e409c4b
Show file tree
Hide file tree
Showing 12 changed files with 412 additions and 1,539 deletions.
120 changes: 58 additions & 62 deletions scripts/evansPolanyi.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,46 +2,39 @@
# encoding: utf-8

"""
This script will generate an Evans-Polanyi plot for a single kinetics
depository.
This script will generate an Evans-Polanyi plot for a single kinetics family.
"""

import math
import numpy
import pylab

from rmgpy.species import Species
from rmgpy.reaction import Reaction
from rmgpy.kinetics import Arrhenius

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

def generateThermoData(species, database):
def generate_thermo_data(species, database):
"""
Return the thermodynamics data for a given `species` as generated from
the specified `database`.
"""
thermoData = [database.thermo.getThermoData(molecule) for molecule in species.molecule]
thermoData.sort(key=lambda x: x.getEnthalpy(298))
return thermoData[0]
thermo_data = [database.thermo.get_thermo_data(species)]
thermo_data.sort(key=lambda x: x.get_enthalpy(298))
return thermo_data[0]

def getReactionForEntry(entry, database):
def get_reaction_for_entry(entry, database):
"""
Return a Reaction object for a given entry that uses Species instead of
Molecules (so that we can compute the reaction thermo).
"""
reaction = Reaction(reactants=[], products=[])
for molecule in entry.item.reactants:
molecule.makeHydrogensExplicit()
reactant = Species(molecule=[molecule], label=molecule.toSMILES())
for reactant in entry.item.reactants:
reactant.generate_resonance_structures()
reactant.thermo = generateThermoData(reactant, database)
reactant.thermo = generate_thermo_data(reactant, database)
reaction.reactants.append(reactant)
for molecule in entry.item.products:
molecule.makeHydrogensExplicit()
product = Species(molecule=[molecule], label=molecule.toSMILES())
for product in entry.item.products:
product.generate_resonance_structures()
product.thermo = generateThermoData(product, database)
product.thermo = generate_thermo_data(product, database)
reaction.products.append(product)

reaction.kinetics = entry.data
Expand All @@ -51,7 +44,7 @@ def getReactionForEntry(entry, database):

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

def fitEvansPolanyi(dHrxn, Ea):
def fit_evans_polanyi(dHrxn, Ea):
"""
Generate best-fit Evans-Polanyi parameters for the given set of enthalpy
of reaction and activation energy data. This can be done using a simple
Expand All @@ -66,7 +59,7 @@ def fitEvansPolanyi(dHrxn, Ea):
A[index,0] = dHrxn[index]
A[index,1] = 1
b[index] = Ea[index]
x, residues, rank, s = numpy.linalg.lstsq(A, b)
x, residues, rank, s = numpy.linalg.lstsq(A, b, rcond=None)

# In the above fit, we probably included reactions with very negative or
# very positive dHrxn, for which Ea = 0 and Ea = dHrxn, respectively
Expand All @@ -77,16 +70,18 @@ def fitEvansPolanyi(dHrxn, Ea):
# excluding the very negative and very positive dHrxn points from our
# dataset and refitting using only the points in the linear region
# This presumes that we will impose 0 <= Ea <= dHrxn where needed
Hmin = -x[1] / x[0]; Hmax = x[1] / (1 - x[0])
Hmin = -x[1] / x[0]
Hmax = x[1] / (1 - x[0])
A = []; b = []
for index in range(N):
if Hmin <= dHrxn[index] <= Hmax:
A.append([dHrxn[index], 1])
b.append(Ea[index])
A = numpy.array(A, numpy.float64)
b = numpy.array(b, numpy.float64)
x, residues, rank, s = numpy.linalg.lstsq(A, b)
Hmin = -x[1] / x[0]; Hmax = x[1] / (1 - x[0])
x, residues, rank, s = numpy.linalg.lstsq(A, b, rcond=None)
Hmin = -x[1] / x[0]
Hmax = x[1] / (1 - x[0])

# Compute sample standard deviation and standard error
stdev = 0.0; error = 0.0; count = 0
Expand Down Expand Up @@ -117,7 +112,7 @@ def fitEvansPolanyi(dHrxn, Ea):

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

def generateEvansPolanyiPlot(depository, database):
def generate_evans_polanyi_plot(depository, database):
"""
Generate an Evans-Polanyi plot for the entries in the given `depository`
of the loaded `database`.
Expand All @@ -126,29 +121,28 @@ def generateEvansPolanyiPlot(depository, database):
fig = pylab.figure(figsize=(6,5))

Ea = []; dHrxn = []; reactions = []
entries = depository.entries.values()
entries = list(depository.entries.values())

for entry in entries:
if isinstance(entry.data, Arrhenius):

reaction = getReactionForEntry(entry, database)
reaction = get_reaction_for_entry(entry, database)

reactions.append(reaction)
Ea.append(reaction.kinetics.Ea.value / 1000.)
dHrxn.append(reaction.getEnthalpyOfReaction(298) / 1000.)
dHrxn.append(reaction.get_enthalpy_of_reaction(298) / 1000.)

xEP, stdevEP, dHrxnEP, EaEP = fitEvansPolanyi(dHrxn, Ea)
alpha, E0 = xEP

print
print 'Summary'
print '======='
print 'Parameters:'
print ' Evans-Polanyi = %s' % (xEP)
print 'Sample standard deviation:'
print ' Evans-Polanyi = %g kJ/mol' % (stdevEP)
print '95% confidence limit:'
print ' Evans-Polanyi = %g kJ/mol' % (1.96 * stdevEP)
xEP, stdevEP, dHrxnEP, EaEP = fit_evans_polanyi(dHrxn, Ea)

print()
print('Summary')
print('=======')
print('Parameters:')
print(' Evans-Polanyi = %s' % (xEP))
print('Sample standard deviation:')
print(' Evans-Polanyi = %g kJ/mol' % (stdevEP))
print('95% confidence limit:')
print(' Evans-Polanyi = %g kJ/mol' % (1.96 * stdevEP))

pylab.plot(dHrxn, Ea, color='#B0B0B0', linestyle='None', marker='o', picker=5)
pylab.plot(dHrxnEP, EaEP, color='red', linestyle='solid', linewidth=2)
Expand All @@ -163,24 +157,24 @@ def generateEvansPolanyiPlot(depository, database):

fig.subplots_adjust(left=0.14, bottom=0.1, top=0.95, right=0.95, wspace=0.20, hspace=0.20)

def onpick(event):
print 'Pick'
def on_pick(event):
print('Pick')
thisline = event.artist
dHrxn = thisline.get_xdata()
Ea = thisline.get_ydata()
for ind in event.ind:
entry = entries[ind]
reaction = reactions[ind]
print '%i. %s' % (entry.index, entry.label)
print reaction
print 'dHrxn = %.2f kJ/mol' % (dHrxn[ind])
print 'Ea = %.2f kJ/mol' % (Ea[ind])
print('%i. %s' % (entry.index, entry.label))
print(reaction)
print('dHrxn = %.2f kJ/mol' % (dHrxn[ind]))
print('Ea = %.2f kJ/mol' % (Ea[ind]))
for reactant in reaction.reactants:
print '%24s H298 = %g kcal/mol' % (reactant, reactant.thermo.getEnthalpy(300) / 4184.)
print('%24s H298 = %g kcal/mol' % (reactant, reactant.thermo.get_enthalpy(300) / 4184.))
for product in reaction.products:
print '%24s H298 = %g kcal/mol' % (product, product.thermo.getEnthalpy(300) / 4184.)
print('%24s H298 = %g kcal/mol' % (product, product.thermo.get_enthalpy(300) / 4184.))

connection_id = fig.canvas.mpl_connect('pick_event', onpick)
fig.canvas.mpl_connect('pick_event', on_pick)

pylab.show()

Expand All @@ -191,24 +185,26 @@ def onpick(event):
import argparse

parser = argparse.ArgumentParser()
parser.add_argument('depository', metavar='<depository>', type=str, nargs=1, help='the depository to use')
parser.add_argument('kinetics_family', metavar='<family>', type=str, nargs=1, help='the family to use')
parser.add_argument('kinetics_depository', metavar='<kinetics_depository>', type=str, nargs='+',
help='the kineticsDepository to use, e.g., training, NIST')

args = parser.parse_args()

print 'Loading RMG database...'
print('Loading RMG database...')
from rmgpy.data.rmg import RMGDatabase
from rmgpy import settings
database = RMGDatabase()
database.load(settings['database.directory'], kineticsFamilies='all', kineticsDepositories='all')

try:
depositories = [database.kinetics.depository[label] for label in args.depository]
except KeyError:
print e
print 'The specified depository "{0}" was invalid.'.format(label)
quit()

print 'Generating Evans-Polanyi data...'
for depository in depositories:
generateEvansPolanyiPlot(depository, database)

database.load(settings['database.directory'], kinetics_families=args.kinetics_family,
kinetics_depositories=args.kinetics_depository)

family = database.kinetics.families.get(args.kinetics_family[0])

print('Generating Evans-Polanyi data...')
for depository in family.depositories:
try:
generate_evans_polanyi_plot(depository, database)
except KeyError:
print(e)
print('The specified depository "{0}" was invalid.'.format(depository))
quit()
54 changes: 28 additions & 26 deletions scripts/exportKineticsLibraryToChemkin.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,67 +16,69 @@
"""
import argparse
import os

from rmgpy import settings

from rmgpy.data.rmg import RMGDatabase
from rmgpy.data.kinetics.library import LibraryReaction
from rmgpy.chemkin import saveChemkinFile, saveSpeciesDictionary
from rmgpy.chemkin import save_chemkin_file, save_species_dictionary
from rmgpy.rmg.model import Species
from rmgpy import settings


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

if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument('library', metavar='LIBRARYNAME', type=str, nargs=1,
help='the name of the kinetic library to be exported')
args = parser.parse_args()
libraryName = args.library[0]
library_name = args.library[0]

print 'Loading RMG-Py database...'
print('Loading RMG-Py database...')
database = RMGDatabase()
database.load(settings['database.directory'], kineticsFamilies='all', kineticsDepositories='all')
database.load(settings['database.directory'], kinetics_families='all', kinetics_depositories='all')


print 'Loading {0} library'.format(libraryName)
kineticLibrary = database.kinetics.libraries[libraryName]
print('Loading {0} library'.format(library_name))
kinetic_library = database.kinetics.libraries[library_name]

reactionList = []
for index, entry in kineticLibrary.entries.iteritems():
reaction_list = []
for index, entry in kinetic_library.entries.items():
reaction = entry.item
reaction.kinetics = entry.data
library_reaction = LibraryReaction(index=reaction.index,
library_reaction = LibraryReaction(index=reaction.index,
reactants=reaction.reactants,
products=reaction.products,
reversible=reaction.reversible,
kinetics=reaction.kinetics,
library=libraryName
library=library_name
)
reactionList.append(library_reaction)
reaction_list.append(library_reaction)

speciesList = []
species_list = []
index = 0
speciesDict = kineticLibrary.getSpecies(os.path.join(settings['database.directory'],'kinetics', 'libraries', libraryName, 'dictionary.txt'))
for spec in speciesDict.values():
species_dict = kinetic_library.get_species(os.path.join(settings['database.directory'],'kinetics', 'libraries', library_name, 'dictionary.txt'))
for spec in list(species_dict.values()):
index = index + 1
species = Species(molecule = spec.molecule)
species.getThermoData()
species.get_thermo_data()
species.index = index
speciesList.append(species)
species_list.append(species)

for reaction in reactionList:
for reaction in reaction_list:
for reactant in reaction.reactants:
for spec in speciesList:
if reactant.isIsomorphic(spec):
for spec in species_list:
if reactant.is_isomorphic(spec):
reactant.index = spec.index
spec.label = reactant.label
for product in reaction.products:
for spec in speciesList:
if product.isIsomorphic(spec):
for spec in species_list:
if product.is_isomorphic(spec):
product.index = spec.index
spec.label = product.label


print 'Saving the chem.inp and species_dictionary.txt to local directory'
saveChemkinFile('chem.inp', speciesList, reactionList)
saveSpeciesDictionary('species_dictionary.txt', speciesList)
print('Saving the chem.inp and species_dictionary.txt to local directory')
save_chemkin_file('chem.inp', species_list, reaction_list)
save_species_dictionary('species_dictionary.txt', species_list)


51 changes: 0 additions & 51 deletions scripts/exportOldDatabase.py
View file

This file was deleted.

Loading

0 comments on commit e409c4b

Please sign in to comment.