diff --git a/scripts/evansPolanyi.py b/scripts/evansPolanyi.py index c8a25ba2b2..2015608a17 100644 --- a/scripts/evansPolanyi.py +++ b/scripts/evansPolanyi.py @@ -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 @@ -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 @@ -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 @@ -77,7 +70,8 @@ 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: @@ -85,8 +79,9 @@ def fitEvansPolanyi(dHrxn, Ea): 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 @@ -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`. @@ -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) @@ -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() @@ -191,24 +185,26 @@ def onpick(event): import argparse parser = argparse.ArgumentParser() - parser.add_argument('depository', metavar='', type=str, nargs=1, help='the depository to use') + parser.add_argument('kinetics_family', metavar='', type=str, nargs=1, help='the family to use') + parser.add_argument('kinetics_depository', metavar='', 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) - \ No newline at end of file + 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() diff --git a/scripts/exportKineticsLibraryToChemkin.py b/scripts/exportKineticsLibraryToChemkin.py index 87540cc382..74dd5934e0 100644 --- a/scripts/exportKineticsLibraryToChemkin.py +++ b/scripts/exportKineticsLibraryToChemkin.py @@ -16,12 +16,14 @@ """ 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__': @@ -29,54 +31,54 @@ 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) diff --git a/scripts/exportOldDatabase.py b/scripts/exportOldDatabase.py deleted file mode 100644 index ee29be0ba9..0000000000 --- a/scripts/exportOldDatabase.py +++ /dev/null @@ -1,51 +0,0 @@ -#!/usr/bin/env python -# encoding: utf-8 - -""" -This script exports the database to the old RMG-Java format. The script -requires two command-line arguments: the path to the database to import, and -the path to save the old RMG-Java database to. -""" - -import os -import argparse -from rmgpy.data.rmg import RMGDatabase -from rmgpy import settings - - -############################################################################### - - -def export(input, output, database=None): - - print 'Loading the new RMG-Py database...' - if not database: - database = RMGDatabase() - database.load(input, kineticsFamilies='all', kineticsDepositories='all') - - print 'Constructing additional rate rules from kinetics depository...' - for family in database.kinetics.families.values(): - family.addKineticsRulesFromTrainingSet(thermoDatabase=database.thermo) - - print "Deleting thermo library entries with atoms RMG-Java can't understand..." - database.thermo.pruneHeteroatoms(allowed=['C','H','O','S']) - print 'Saving old RMG-Java database...' - database.saveOld(output) - print "Done!" - - -############################################################################### - - - -if __name__ == '__main__': - - parser = argparse.ArgumentParser() - parser.add_argument('outputPath', metavar='OUTPUT', type=str, nargs=1, - help='the outputPath for the RMG-Java database directory') - - args = parser.parse_args() - outputPath = args.outputPath[0] - - inputPath = settings['database.directory'] - export(inputPath, outputPath) diff --git a/scripts/fitPolycyclicThermoGroupsFromThermoLibrary.ipynb b/scripts/fitPolycyclicThermoGroupsFromThermoLibrary.ipynb index 1095df53ff..92bde2a02f 100644 --- a/scripts/fitPolycyclicThermoGroupsFromThermoLibrary.ipynb +++ b/scripts/fitPolycyclicThermoGroupsFromThermoLibrary.ipynb @@ -18,35 +18,27 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "# Fill in the list of thermo libraries to be used for fitting polycyclic thermo groups\n", - "thermoLibraries = ['vinylCPD_H','C3','C10H11','Fulvene_H','naphthalene_H']\n" + "thermo_libraries = ['vinylCPD_H','C3','C10H11','Fulvene_H','naphthalene_H']\n" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ - "import os\n", - "import re\n", "import copy\n", - "import numpy\n", - "from IPython.display import Image, display\n", + "from IPython.display import display\n", "from rmgpy.data.thermo import *\n", "from rmgpy.data.base import Entry\n", "from rmgpy.data.rmg import RMGDatabase\n", "from rmgpy import settings\n", "from rmgpy.species import Species\n", - "from rmgpy.molecule import Molecule\n", - "from rmgpy.cantherm.output import prettify" + "from arkane.output import prettify" ] }, { @@ -59,50 +51,48 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ - "def extractPolycyclicGroups(molecule):\n", + "def extract_polycyclic_groups(molecule):\n", " \"\"\"\n", " Extract polycyclic functional groups from a real molecule\n", " \"\"\"\n", " struct = molecule.copy(deep=True)\n", " # Saturate the structure if it is a radical\n", - " if struct.isRadical():\n", + " if struct.is_radical():\n", " struct.saturate()\n", - " struct.deleteHydrogens()\n", + " struct.delete_hydrogens()\n", " \n", - " polyRings = struct.getPolycyclicRings()\n", - " groups = [convertCycleToGroup(ring) for ring in polyRings]\n", + " poly_rings = struct.get_polycycles()\n", + " groups = [convert_cycle_to_group(ring) for ring in poly_rings]\n", " \n", " return groups\n", " \n", - "def convertCycleToGroup(cycle):\n", + "def convert_cycle_to_group(cycle):\n", " \"\"\"\n", " This function converts a list of atoms in a cycle to a functional Group object\n", " \"\"\"\n", " from rmgpy.molecule.group import GroupAtom, GroupBond, Group\n", " \n", " # Create GroupAtom object for each atom in the cycle, label the first one in the cycle with a *\n", - " groupAtoms = {}\n", + " group_atoms = {}\n", " bonds = []\n", " for atom in cycle:\n", - " groupAtoms[atom] = GroupAtom(atomType=[atom.atomType],\n", - " radicalElectrons=[0],\n", + " group_atoms[atom] = GroupAtom(atomType=[atom.atomtype],\n", + " radical_electrons=[0],\n", " label='*' if cycle.index(atom)==0 else '')\n", " \n", - " group = Group(atoms=groupAtoms.values()) \n", + " group = Group(atoms=group_atoms.values()) \n", " \n", " # Create GroupBond for each bond between atoms in the cycle, but not outside of the cycle\n", " for atom in cycle:\n", - " for bondedAtom, bond in atom.edges.iteritems():\n", - " if bondedAtom in cycle:\n", + " for bonded_atom, bond in atom.edges.iteritems():\n", + " if bonded_atom in cycle:\n", " # create a group bond with the same bond order as in the original molecule,\n", " # if it hasn't already been created\n", - " if not group.hasBond(groupAtoms[atom],groupAtoms[bondedAtom]):\n", - " group.addBond(GroupBond(groupAtoms[atom],groupAtoms[bondedAtom],order=[bond.order]))\n", + " if not group.has_bond(group_atoms[atom],group_atoms[bonded_atom]):\n", + " group.add_bond(GroupBond(group_atoms[atom],group_atoms[bonded_atom],order=[bond.order]))\n", " else:\n", " pass\n", " \n", @@ -121,23 +111,21 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ - "def displayThermo(thermoData):\n", - " print 'H298 = {0} kcal/mol'.format(thermoData.H298.value_si/4184)\n", - " print 'S298 = {0} cal/mol*K'.format(thermoData.S298.value_si/4.184)\n", - "def compareThermoData(thermoData1, thermoData2):\n", - " delH = thermoData1.H298.value_si - thermoData2.H298.value_si\n", - " print 'Difference in H298 = {0} kcal/mol'.format(delH/4184)\n", - " delS = thermoData1.S298.value_si - thermoData2.S298.value_si\n", - " print 'Difference S298 = {0} cal/mol*K'.format(delS/4.184)\n", - " #Tdata = [300,500,1000,2000]\n", - " #for T in Tdata:\n", - " # delCp = thermoData1.getHeatCapacity(T) - thermoData2.getHeatCapacity(T)\n", - " # print 'Difference in Cp at {0} = {1} cal/mol*K'.format(T, delCp/4.184)" + "def display_thermo(thermoData):\n", + " print('H298 = {0} kcal/mol'.format(thermoData.H298.value_si/4184))\n", + " print('S298 = {0} cal/mol*K'.format(thermoData.S298.value_si/4.184))\n", + "def compare_thermo_data(thermo_data1, thermo_data2):\n", + " del_H = thermo_data1.H298.value_si - thermo_data2.H298.value_si\n", + " print('Difference in H298 = {0} kcal/mol'.format(del_H/4184))\n", + " del_S = thermo_data1.S298.value_si - thermo_data2.S298.value_si\n", + " print('Difference S298 = {0} cal/mol*K'.format(del_S/4.184))\n", + " #T_data = [300,500,1000,2000]\n", + " #for T in T_data:\n", + " # del_cp = thermo_data1.get_heat_capacity(T) - thermo_data2.get_heat_capacity(T)\n", + " # print('Difference in cp at {0} = {1} cal/mol*K'.format(T, del_cp/4.184))" ] }, { @@ -150,16 +138,14 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "database = RMGDatabase()\n", - "database.load(settings['database.directory'], thermoLibraries = thermoLibraries, kineticsFamilies='none', kineticsDepositories='none', reactionLibraries=[])\n", + "database.load(settings['database.directory'], thermo_libraries = thermo_libraries, kinetics_families='none', kinetics_depositories='none', reaction_libraries=[])\n", "\n", - "thermoDatabase = database.thermo\n", - "thermoDatabase0 = copy.deepcopy(database.thermo)" + "thermo_database = database.thermo\n", + "thermo_database0 = copy.deepcopy(database.thermo)" ] }, { @@ -174,93 +160,93 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ - "fittingDictionary={}\n", - "for libraryName in thermoLibraries:\n", - " thermoLibrary = database.thermo.libraries[libraryName]\n", - " for label, entry in thermoLibrary.entries.iteritems():\n", + "fitting_dictionary={}\n", + "for library_name in thermo_libraries:\n", + " thermo_library = database.thermo.libraries[library_name]\n", + " for label, entry in iter(thermo_library.entries.items()):\n", " molecule = entry.item\n", - " libraryThermoData = entry.data\n", - " if molecule.getAllPolycyclicVertices():\n", - " print label\n", + " library_thermo_data = entry.data\n", + " if molecule.get_all_polycyclic_vertices():\n", + " print(label)\n", " species = Species(molecule=[molecule])\n", " species.generate_resonance_structures() \n", - " print 'Species has {0} resonance structures'.format(len(species.molecule))\n", - " estimatedThermos = [thermoDatabase.estimateThermoViaGroupAdditivity(molecule) for molecule in species.molecule]\n", - " for estimatedThermo in estimatedThermos:\n", - " index = estimatedThermos.index(estimatedThermo)\n", - " ringGroups, polycyclicGroups = database.thermo.getRingGroupsFromComments(estimatedThermo)\n", + " print('Species has {0} resonance structures'.format(len(species.molecule)))\n", + " estimated_thermos = [thermo_database.estimate_thermo_via_group_additivity(molecule) \n", + " for molecule in species.molecule]\n", + " for estimated_thermo in estimated_thermos:\n", + " index = estimated_thermos.index(estimated_thermo)\n", + " ring_groups, polycyclic_groups = database.thermo.get_ring_groups_from_comments(estimated_thermo)\n", " \n", - " if len(polycyclicGroups) == 0:\n", - " raise Exception('Species {0} detected as polycyclic but estimated thermo contained no polycyclic groups: \\\n", - " you need to create a new polycyclic group'.format(label))\n", + " if len(polycyclic_groups) == 0:\n", + " raise Exception('Species {0} detected as polycyclic but estimated thermo contained no \\\n", + " polycyclic groups: you need to create a new polycyclic group'.format(label))\n", "\n", - " elif len(polycyclicGroups) == 1:\n", - " polycyclicGroup = polycyclicGroups[0]\n", - " print 'Molecule {0} has a single polycyclic group match in thermo estimate.'.format(species.molecule[index].toSMILES())\n", + " elif len(polycyclic_groups) == 1:\n", + " polycyclic_group = polycyclic_groups[0]\n", + " print('Molecule {0} has a single polycyclic group match in thermo estimate.'.format(\n", + " species.molecule[index].to_smiles()))\n", " # Draw the molecule in ipython notebook\n", " display(species.molecule[index])\n", - " print 'Molecule SMILES: {0}'.format(species.molecule[index].toSMILES())\n", - " print 'Estimated thermo data:'\n", - " print prettify(repr(estimatedThermo))\n", - " displayThermo(estimatedThermo)\n", + " print('Molecule SMILES: {0}'.format(species.molecule[index].to_smiles()))\n", + " print('Estimated thermo data:')\n", + " print(prettify(repr(estimated_thermo)))\n", + " display_thermo(estimated_thermo)\n", "\n", - " withoutPolycyclicGroupThermo = removeThermoData(copy.deepcopy(estimatedThermo), polycyclicGroup.data)\n", - " newPolycyclicGroupThermo = removeThermoData(copy.deepcopy(libraryThermoData), withoutPolycyclicGroupThermo)\n", + " without_polycyclic_group_thermo = remove_thermo_data(copy.deepcopy(estimated_thermo), \n", + " polycyclic_group.data)\n", + " new_polycyclic_group_thermo = remove_thermo_data(copy.deepcopy(library_thermo_data), \n", + " without_polycyclic_group_thermo)\n", "\n", "\n", " # Check to make sure that the polycyclic group is not generic\n", " # If it is, create a new polycyclicGroup as the child\n", - " if polycyclicGroup.label == 'PolycyclicRing':\n", - " groups = extractPolycyclicGroups(species.molecule[index])\n", - " print groups[0].toAdjacencyList()\n", + " if polycyclic_group.label == 'PolycyclicRing':\n", + " groups = extract_polycyclic_groups(species.molecule[index])\n", + " print(groups[0].to_adjacency_list())\n", " assert len(groups) == 1\n", " # Create a new entry in the polycyclic groups with the same name as the thermo library entry\n", " # Make sure it does not already exist\n", - " entryLabel = label\n", + " entry_label = label\n", " counter = 0\n", - " while entryLabel in thermoDatabase.groups['polycyclic'].entries:\n", + " while entry_label in thermoDatabase.groups['polycyclic'].entries:\n", " counter += 1\n", - " entryLabel = entryLabel.split('-')[0]\n", - " entryLabel += '-{0}'.format(counter)\n", - "\n", + " entry_label = entry_label.split('-')[0]\n", + " entry_label += '-{0}'.format(counter)\n", "\n", - " print 'Polycyclic group was found to be generic \"PolycyclicRing\". Creating new entry: {0}.'.format(entryLabel)\n", - " thermoDatabase.groups['polycyclic'].entries[entryLabel] = Entry(index = len(thermoDatabase.groups['polycyclic'].entries)+1,\n", - " label = entryLabel,\n", - " item = groups[0],\n", - " data = polycyclicGroup.data, # Use dummy thermo here so other estimates can find this group\n", - " parent = polycyclicGroup,\n", - " )\n", "\n", - " # Set the new entry as the polycyclicGroup and make it a child of the generic group\n", - " polycyclicGroup = thermoDatabase.groups['polycyclic'].entries[entryLabel] \n", - " thermoDatabase.groups['polycyclic'].entries['PolycyclicRing'].children.append(polycyclicGroup)\n", + " print('Polycyclic group was found to be generic \"PolycyclicRing\". \\\n", + " Creating new entry: {0}.'.format(entry_label))\n", + " thermo_database.groups['polycyclic'].entries[entry_label] = Entry(index = len(\n", + " thermo_database.groups['polycyclic'].entries)+1,\n", + " label = entry_label,\n", + " item = groups[0],\n", + " data = polycyclic_group.data, # Use dummy thermo here so other estimates can find this group\n", + " parent = polycyclic_group,\n", + " )\n", "\n", + " # Set the new entry as the polycyclic_group and make it a child of the generic group\n", + " polycyclic_group = thermo_database.groups['polycyclic'].entries[entryLabel] \n", + " thermo_database.groups['polycyclic'].entries['PolycyclicRing'].children.append(\n", + " polycyclic_group)\n", "\n", " else:\n", - " print 'Matched polycyclic group \"{0}\"'.format(polycyclicGroup.label)\n", - "\n", - "\n", - "\n", - "\n", + " print('Matched polycyclic group \"{0}\"'.format(polycyclic_group.label))\n", "\n", " # Add the new group value to the fitting dictionary\n", - " if polycyclicGroup not in fittingDictionary:\n", + " if polycyclic_group not in fitting_dictionary:\n", " # Add a tuple containing fitted group data, the original library entry, and thermo library\n", - " fittingDictionary[polycyclicGroup]=[(newPolycyclicGroupThermo, entry, thermoLibrary)]\n", + " fitting_dictionary[polycyclic_group]=[(new_polycyclic_group_thermo, entry, thermo_library)]\n", " else:\n", - " fittingDictionary[polycyclicGroup].append((newPolycyclicGroupThermo, entry, thermoLibrary))\n", + " fitting_dictionary[polycyclic_group].append((new_polycyclic_group_thermo, entry, \n", + " thermo_library))\n", "\n", - " elif len(polycyclicGroups) > 1:\n", - " print 'Species {0} has matched multiple polycyclic groups. \\\n", - " This cannot be fitted with a single molecule\\'s thermo data.'.format(label)\n", + " elif len(polycyclic_groups) > 1:\n", + " print('Species {0} has matched multiple polycyclic groups. This cannot be fitted with a single molecule\\'s thermo data.'.format(label))\n", " raise Exception\n", - " print '====================='" + " print('=====================')" ] }, { @@ -273,40 +259,38 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ - "for polycyclicGroup, fittingGroups in fittingDictionary.iteritems():\n", - " print 'Original thermo data for polycyclic group: {0}'.format(polycyclicGroup.label)\n", - " if polycyclicGroup.data:\n", - " print prettify(repr(polycyclicGroup.data))\n", + "for polycyclic_group, fitting_groups in iter(fitting_dictionary.items()):\n", + " print('Original thermo data for polycyclic group: {0}'.format(polycyclic_group.label))\n", + " if polycyclic_group.data:\n", + " print(prettify(repr(polycyclic_group.data)))\n", " else:\n", - " print 'No data found. Was created as a new group.'\n", + " print('No data found. Was created as a new group.')\n", " \n", - " thermoDataset = [fitTuple[0] for fitTuple in fittingGroups]\n", - " labels = [fitTuple[1].label for fitTuple in fittingGroups]\n", - " libraryLabels = [fitTuple[2].name for fitTuple in fittingGroups]\n", + " thermo_dataset = [fit_tuple[0] for fit_tuple in fitting_groups]\n", + " labels = [fit_tuple[1].label for fit_tuple in fitting_groups]\n", + " library_labels = [fit_tuple[2].name for fit_tuple in fitting_groups]\n", " # Average the new group values to fit the original polycyclic group\n", - " fittedGroupData = averageThermoData([fitTuple[0] for fitTuple in fittingGroups])\n", - " #print fittedGroupData\n", - " #print fittingGroups\n", - " polycyclicGroup.data = fittedGroupData\n", - " polycyclicGroup.shortDesc = \"Fitted from thermo library values\"\n", + " fitted_group_data = average_thermo_data([fit_tuple[0] for fit_tuple in fitting_groups])\n", + " #print fitted_group_data\n", + " #print fitting_groups\n", + " polycyclic_group.data = fitted_group_data\n", + " polycyclic_group.short_desc = \"Fitted from thermo library values\"\n", " \n", " comment = ''\n", " for i in range(len(labels)):\n", - " comment += \"Fitted from species {0} from {1} library.\\n\".format(labels[i],libraryLabels[i])\n", - " polycyclicGroup.longDesc = comment.strip()\n", + " comment += \"Fitted from species {0} from {1} library.\\n\".format(labels[i],library_labels[i])\n", + " polycyclic_group.long_desc = comment.strip()\n", " \n", - " print 'Fitted thermo data for polycyclic group: {0}'.format(polycyclicGroup.label)\n", - " print prettify(repr(polycyclicGroup.data))\n", - " print polycyclicGroup.longDesc\n", - " print '===================================================================='\n", + " print('Fitted thermo data for polycyclic group: {0}'.format(polycyclic_group.label))\n", + " print(prettify(repr(polycyclic_group.data)))\n", + " print(polycyclic_group.long_desc)\n", + " print('====================================================================')\n", " # At this point, save and overwrite the entire polycyclic thermo library\n", "\n", - "thermoDatabase.groups['polycyclic'].save('new_polycyclic.py')" + "thermo_database.groups['polycyclic'].save('new_polycyclic.py')" ] }, { @@ -319,96 +303,112 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Test that the new group additivity values can estimate the old library ones.\n", - "\n", - "for libraryName in thermoLibraries:\n", - " thermoLibrary = database.thermo.libraries[libraryName]\n", - " for label, entry in thermoLibrary.entries.iteritems():\n", + "for library_name in thermo_libraries:\n", + " thermo_library = database.thermo.libraries[library_name]\n", + " for label, entry in iter(thermo_library.entries.items()):\n", " molecule = entry.item\n", - " libraryThermoData = entry.data\n", + " library_thermo_data = entry.data\n", "\n", - " if molecule.getAllPolycyclicVertices():\n", + " if molecule.get_all_polycyclic_vertices():\n", " species = Species(molecule=[molecule])\n", " species.generate_resonance_structures()\n", - " print label\n", + " print(label)\n", " display (species.molecule[0])\n", - " print 'Species has {0} resonance structures(s)'.format(len(species.molecule))\n", - " thermoDatabase.findCp0andCpInf(species, libraryThermoData)\n", + " print('Species has {0} resonance structures(s)'.format(len(species.molecule)))\n", + " find_cp0_and_cpinf(species, library_thermo_data)\n", " \n", - " estimatedThermo = thermoDatabase.getThermoDataFromGroups(species)\n", - " originalEstimatedThermo = thermoDatabase0.getThermoDataFromGroups(species)\n", - " if libraryThermoData.isIdenticalTo(estimatedThermo):\n", - " print 'Library thermo data for species {0} matches the estimate from group additivity.\\n'.format(label)\n", + " estimated_thermo = thermo_database.get_thermo_data_from_groups(species)\n", + " original_estimated_thermo = thermo_database0.get_thermo_data_from_groups(species)\n", + " if library_thermo_data.is_identical_to(estimated_thermo):\n", + " print('Library thermo data for species {0} matches the estimate from group additivity.\\n'.format(\n", + " label))\n", " \n", - " print 'Library thermo data:' \n", - " displayThermo(libraryThermoData)\n", - " print '' \n", + " print('Library thermo data:') \n", + " display_thermo(library_thermo_data)\n", + " print('') \n", " \n", - " print 'Original estimated thermo data:' \n", - " ringGroups, polycyclicGroups = database.thermo.getRingGroupsFromComments(estimatedThermo)\n", - " print 'Polycyclic groups: {0}'.format(' '.join([grp.label for grp in polycyclicGroups]))\n", - " displayThermo(originalEstimatedThermo)\n", - " print ''\n", + " print('Original estimated thermo data:') \n", + " ring_groups, polycyclic_groups = database.thermo.get_ring_groups_from_comments(estimated_thermo)\n", + " print('Polycyclic groups: {0}'.format(' '.join([grp.label for grp in polycyclic_groups])))\n", + " display_thermo(original_estimated_thermo)\n", + " print('')\n", " \n", - " print 'Comparison of library thermo with original estimated thermo'\n", - " compareThermoData(libraryThermoData,originalEstimatedThermo)\n", - " print ''\n", + " print('Comparison of library thermo with original estimated thermo')\n", + " compare_thermo_data(library_thermo_data, original_estimated_thermo)\n", + " print('')\n", " else:\n", - " print 'Library thermo data for species {0} does not match the estimate from group additivity\\n'.format(label)\n", + " print('Library thermo data for species {0} does not match the estimate from group \\\n", + " additivity\\n'.format(label))\n", " \n", - " print 'Library thermo data:' \n", - " displayThermo(libraryThermoData)\n", - " print '' \n", + " print('Library thermo data:')\n", + " display_thermo(library_thermo_data)\n", + " print('') \n", " \n", - " print 'Estimated thermo data:' \n", + " print('Estimated thermo data:') \n", " \n", - " ringGroups, polycyclicGroups = database.thermo.getRingGroupsFromComments(estimatedThermo)\n", - " print 'Polycyclic groups: {0}'.format(' '.join([grp.label for grp in polycyclicGroups]))\n", - " displayThermo(estimatedThermo)\n", - " print ''\n", + " ring_groups, polycyclic_groups = database.thermo.get_ring_groups_from_comments(estimated_thermo)\n", + " print('Polycyclic groups: {0}'.format(' '.join([grp.label for grp in polycyclic_groups])))\n", + " display_thermo(estimated_thermo)\n", + " print('')\n", " \n", - " print 'Comparison of library thermo with estimated thermo'\n", - " compareThermoData(libraryThermoData,estimatedThermo)\n", - " print ''\n", + " print('Comparison of library thermo with estimated thermo')\n", + " compare_thermo_data(library_thermo_data, estimated_thermo)\n", + " print('')\n", " \n", - " print 'Original estimated thermo data:' \n", + " print('Original estimated thermo data:') \n", " \n", - " ringGroups, polycyclicGroups = database.thermo.getRingGroupsFromComments(originalEstimatedThermo)\n", - " print 'Polycyclic groups: {0}'.format(' '.join([grp.label for grp in polycyclicGroups]))\n", - " displayThermo(originalEstimatedThermo)\n", - " print ''\n", + " ring_groups, polycyclic_groups = database.thermo.get_ring_groups_from_comments(\n", + " original_estimated_thermo)\n", + " print('Polycyclic groups: {0}'.format(' '.join([grp.label for grp in polycyclic_groups])))\n", + " display_thermo(original_estimated_thermo)\n", + " print('')\n", " \n", - " print 'Comparison of library thermo with original estimated thermo'\n", - " compareThermoData(libraryThermoData,originalEstimatedThermo)\n", - " print ''\n", - " print '======================='" + " print('Comparison of library thermo with original estimated thermo')\n", + " compare_thermo_data(library_thermo_data, original_estimated_thermo)\n", + " print('')\n", + " print('=======================')" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python 2", + "display_name": "Python 3", "language": "python", - "name": "python2" + "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", - "version": 2 + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.11" + "pygments_lexer": "ipython3", + "version": "3.7.4" + }, + "pycharm": { + "stem_cell": { + "cell_type": "raw", + "source": [], + "metadata": { + "collapsed": false + } + } } }, "nbformat": 4, - "nbformat_minor": 0 -} + "nbformat_minor": 1 +} \ No newline at end of file diff --git a/scripts/generateFilterArrheniusFits.ipynb b/scripts/generateFilterArrheniusFits.ipynb index 1ebc4fb584..2fc94a3141 100644 --- a/scripts/generateFilterArrheniusFits.ipynb +++ b/scripts/generateFilterArrheniusFits.ipynb @@ -23,8 +23,6 @@ "metadata": {}, "outputs": [], "source": [ - "from __future__ import division, print_function\n", - "\n", "import os\n", "import operator\n", "\n", @@ -33,7 +31,7 @@ "import numpy\n", "\n", "from rmgpy import settings\n", - "from rmgpy.data.kinetics.database import FilterLimitFits\n", + "from rmgpy.data.kinetics.database import filter_limit_fits\n", "from rmgpy.data.rmg import RMGDatabase\n", "from rmgpy.kinetics.arrhenius import Arrhenius\n", "from rmgpy.thermo.thermoengine import submit" @@ -55,7 +53,7 @@ "database = RMGDatabase()\n", "database.load(\n", " settings['database.directory'], \n", - " thermoLibraries = [\n", + " thermo_libraries = [\n", " 'primaryThermoLibrary',\n", " 'Klippenstein_Glarborg2016',\n", " 'BurkeH2O2',\n", @@ -67,11 +65,11 @@ " 'SABIC_aromatics',\n", " 'vinylCPD_H'\n", " ],\n", - " transportLibraries = [],\n", - " reactionLibraries = [],\n", - " seedMechanisms = [],\n", - " kineticsFamilies = 'all',\n", - " kineticsDepositories = ['training'],\n", + " transport_libraries = [],\n", + " reaction_libraries = [],\n", + " seed_mechanisms = [],\n", + " kinetics_families = 'all',\n", + " kinetics_depositories = ['training'],\n", " depository = False, \n", ")" ] @@ -97,10 +95,10 @@ "outputs": [], "source": [ "# Temperature range to fit Arrhenius\n", - "Tmin = 298.0\n", - "Tmax = 2500.0\n", - "Tcount = 50\n", - "Ts = 1 / numpy.linspace(1 / Tmax, 1 / Tmin, Tcount)" + "T_min = 298.0\n", + "T_max = 2500.0\n", + "T_count = 50\n", + "Ts = 1 / numpy.linspace(1 / T_max, 1 / T_min, T_count)" ] }, { @@ -121,7 +119,7 @@ " print(fam_name)\n", " \n", " fam = database.kinetics.families[fam_name]\n", - " dep = fam.getTrainingDepository()\n", + " dep = fam.get_training_depository()\n", " rxns = []\n", "\n", " # Extract all training reactions for selected family\n", @@ -144,7 +142,7 @@ " k_list.append(rxn.kinetics)\n", " index_list.append(rxn.index)\n", " if len(rxn.products) == molecularity:\n", - " k_list.append(rxn.generateReverseRateCoefficient())\n", + " k_list.append(rxn.generate_reverse_rate_coefficient())\n", " index_list.append(rxn.index)\n", "\n", " # Get max. kinetic rates at each discrete temperature\n", @@ -152,30 +150,30 @@ " k_max_list = []\n", " max_rxn_list = set()\n", " for T in Ts:\n", - " kvals = [k.getRateCoefficient(T) for k in k_list]\n", + " kvals = [k.get_rate_coefficient(T) for k in k_list]\n", " mydict = dict(zip(index_list, kvals))\n", "\n", " # Find key and value of max rate coefficient\n", - " key_max_rate = max(mydict.iteritems(), key=operator.itemgetter(1))[0]\n", + " key_max_rate = max(mydict.items(), key=operator.itemgetter(1))[0]\n", " \n", " max_entry = dep.entries.get(key_max_rate)\n", " max_rxn = max_entry.item\n", " max_rxn_list.add(max_rxn)\n", " \n", - " kval = mydict[key_max_rate]\n", - " k_max_list.append(kval)\n", + " k_val = mydict[key_max_rate]\n", + " k_max_list.append(k_val)\n", " \n", - " if molecularity==2 and max_rxn.check_collision_limit_violation(Tmin, Tmax, 10000.0, 1.0e7):\n", + " if molecularity==2 and max_rxn.check_collision_limit_violation(T_min, T_max, 10000.0, 1.0e7):\n", " display(max_rxn)\n", " print(\"\"\"The collision limit of {0} m^3/(mol*s) at {1} (K) is violated by \n", " training reaction {2} with index {3}.\n", " \n", " The rate of training reaction {2} \n", - " is {4} m^3/(mol*s).\"\"\".format(max_rxn.calculate_coll_limit(T), T, max_entry, key_max_rate, kval))\n", + " is {4} m^3/(mol*s).\"\"\".format(max_rxn.calculate_coll_limit(T), T, max_entry, key_max_rate, k_val))\n", "\n", " units = 's^-1' if molecularity == 1 else 'm^3/(mol*s)'\n", " \n", - " arr = Arrhenius().fitToData(Ts, numpy.array(k_max_list), units)\n", + " arr = Arrhenius().fit_to_data(Ts, numpy.array(k_max_list), units)\n", " \n", " fig = plt.figure()\n", " fig_name = fam_name\n", @@ -259,7 +257,7 @@ "outputs": [], "source": [ "# Generate empty ArrheniusRMGObject\n", - "filter_fits = FilterLimitFits(unimol=dict_unimol, bimol=dict_bimol)" + "filter_fits = filter_limit_fits(unimol=dict_unimol, bimol=dict_bimol)" ] }, { @@ -283,32 +281,32 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 2", + "display_name": "Python 3", "language": "python", - "name": "python2" + "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", - "version": 2 + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.16" + "pygments_lexer": "ipython3", + "version": "3.7.4" }, "pycharm": { "stem_cell": { "cell_type": "raw", + "source": [], "metadata": { "collapsed": false - }, - "source": [] + } } } }, "nbformat": 4, "nbformat_minor": 2 -} +} \ No newline at end of file diff --git a/scripts/importChemkinLibrary.py b/scripts/importChemkinLibrary.py index 6742809798..1e8c349af9 100644 --- a/scripts/importChemkinLibrary.py +++ b/scripts/importChemkinLibrary.py @@ -10,52 +10,48 @@ import argparse import logging import os -from rmgpy.data.thermo import ThermoLibrary -from rmgpy.data.kinetics import KineticsLibrary -from rmgpy.data.base import Entry -from rmgpy.chemkin import loadChemkinFile, getSpeciesIdentifier + from rmgpy import settings - +from rmgpy.chemkin import load_chemkin_file +from rmgpy.data.base import Entry +from rmgpy.data.kinetics import KineticsLibrary +from rmgpy.data.thermo import ThermoLibrary + if __name__ == '__main__': parser = argparse.ArgumentParser() - parser.add_argument('chemkinPath', metavar='CHEMKIN', type=str, nargs=1, + parser.add_argument('chemkin_path', metavar='CHEMKIN', type=str, nargs=1, help='The path of the chemkin file') - parser.add_argument('dictionaryPath', metavar='DICTIONARY', type=str, nargs=1, + parser.add_argument('dictionary_path', metavar='DICTIONARY', type=str, nargs=1, help='The path of the RMG dictionary file') parser.add_argument('name', metavar='NAME', type=str, nargs=1, help='Name of the chemkin library to be saved') args = parser.parse_args() - chemkinPath = args.chemkinPath[0] - dictionaryPath = args.dictionaryPath[0] + chemkin_path = args.chemkin_path[0] + dictionary_path = args.dictionary_path[0] name = args.name[0] - speciesList, reactionList = loadChemkinFile(chemkinPath, dictionaryPath) + species_list, reaction_list = load_chemkin_file(chemkin_path, dictionary_path) - # Make full species identifier the species labels - for species in speciesList: - species.label = getSpeciesIdentifier(species) - species.index = -1 - # load thermo library entries - thermoLibrary = ThermoLibrary(name=name) - for i in range(len(speciesList)): - species = speciesList[i] + thermo_library = ThermoLibrary(name=name) + for i in range(len(species_list)): + species = species_list[i] if species.thermo: - thermoLibrary.loadEntry(index = i + 1, + thermo_library.load_entry(index = i + 1, label = species.label, - molecule = species.molecule[0].toAdjacencyList(), + molecule = species.molecule[0].to_adjacency_list(), thermo = species.thermo, - shortDesc = species.thermo.comment - ) + ) else: - logging.warning('Species {0} did not contain any thermo data and was omitted from the thermo library.'.format(str(species))) + logging.warning("""Species {0} did not contain any thermo data and was omitted from the thermo + library.""".format(str(species))) # load kinetics library entries - kineticsLibrary = KineticsLibrary(name=name) - kineticsLibrary.entries = {} - for i in range(len(reactionList)): - reaction = reactionList[i] + kinetics_library = KineticsLibrary(name=name) + kinetics_library.entries = {} + for i in range(len(reaction_list)): + reaction = reaction_list[i] entry = Entry( index = i+1, label = str(reaction), @@ -63,24 +59,24 @@ data = reaction.kinetics, ) try: - entry.longDesc = 'Originally from reaction library: ' + reaction.library + "\n" + reaction.kinetics.comment - except AttributeError: - entry.longDesc = reaction.kinetics.comment - kineticsLibrary.entries[i+1] = entry + entry.long_desc = 'Originally from reaction library: ' + reaction.library + "\n" + reaction.kinetics.comment + except AttributeError: + entry.long_desc = reaction.kinetics.comment + kinetics_library.entries[i+1] = entry # Mark as duplicates where there are mixed pressure dependent and non-pressure dependent duplicate kinetics # Even though CHEMKIN does not require a duplicate flag, RMG needs it. - # Using flag markDuplicates = True - kineticsLibrary.checkForDuplicates(markDuplicates=True) - kineticsLibrary.convertDuplicatesToMulti() + # Using flag mark_duplicates = True + kinetics_library.check_for_duplicates(mark_duplicates=True) + kinetics_library.convert_duplicates_to_multi() # Save in Py format - databaseDirectory = settings['database.directory'] + database_directory = settings['database.directory'] try: - os.makedirs(os.path.join(databaseDirectory, 'kinetics', 'libraries',name)) + os.makedirs(os.path.join(database_directory, 'kinetics', 'libraries',name)) except: pass - thermoLibrary.save(os.path.join(databaseDirectory, 'thermo' ,'libraries', name + '.py')) - kineticsLibrary.save(os.path.join(databaseDirectory, 'kinetics', 'libraries', name, 'reactions.py')) - kineticsLibrary.saveDictionary(os.path.join(databaseDirectory, 'kinetics', 'libraries', name, 'dictionary.txt')) + thermo_library.save(os.path.join(database_directory, 'thermo' ,'libraries', name + '.py')) + kinetics_library.save(os.path.join(database_directory, 'kinetics', 'libraries', name, 'reactions.py')) + kinetics_library.save_dictionary(os.path.join(database_directory, 'kinetics', 'libraries', name, 'dictionary.txt')) diff --git a/scripts/importJavaKineticsLibrary.py b/scripts/importJavaKineticsLibrary.py deleted file mode 100644 index 9a5ed497bc..0000000000 --- a/scripts/importJavaKineticsLibrary.py +++ /dev/null @@ -1,46 +0,0 @@ -#!/usr/bin/env python -# encoding: utf-8 - -""" -This script imports an individual RMG-Java kinetics library from a local directory and saves the output -kinetics library py file into a path of the user's choosing. This library will be automatically -added to the 'libraryname' folder in the input/kinetics/libraries directory and can -be used directly as an RMG-Py kinetics library. - -usage: -importJavaKineticsLibrary.py [-h] INPUT LIBRARYNAME - -positional arguments: -INPUT the input path of the RMG-Java kinetics library directory -LIBRARYNAME the libraryname for the RMG-Py format kinetics library - -""" - -import argparse -import os -from rmgpy.data.kinetics import KineticsLibrary -from rmgpy import settings - -if __name__ == '__main__': - - parser = argparse.ArgumentParser() - parser.add_argument('inputPath', metavar='INPUT', type=str, nargs=1, - help='the input path of the RMG-Java kinetics library directory') - parser.add_argument('libraryName', metavar='OUTPUT', type=str, nargs=1, - help='the libraryName for the RMG-Py format kinetics library') - - args = parser.parse_args() - inputPath = args.inputPath[0] - libraryName = args.libraryName[0] - - library = KineticsLibrary() - library.loadOld(inputPath) - - try: - os.makedirs(os.path.join(settings['database.directory'], 'kinetics', 'libraries', libraryName)) - except: - pass - - # Save in Py format - library.save(os.path.join(settings['database.directory'], 'kinetics', 'libraries', libraryName, 'reactions.py')) - library.saveDictionary(os.path.join(settings['database.directory'], 'kinetics', 'libraries', libraryName,'dictionary.txt')) diff --git a/scripts/importJavaThermoLibrary.py b/scripts/importJavaThermoLibrary.py deleted file mode 100644 index 37ee2ed6c8..0000000000 --- a/scripts/importJavaThermoLibrary.py +++ /dev/null @@ -1,48 +0,0 @@ -#!/usr/bin/env python -# encoding: utf-8 - -""" -This script imports an individual RMG-Java themo library from a local directory and saves the output -thermo library py file into a path of the user's choosing. This library will be automatically -saved to libraryname.py in the input/thermo/libraries directory and can -be used directly as an RMG-Py thermo library. - -usage: -importJavaThermoLibrary.py [-h] INPUT LIBRARYNAME - -positional arguments: -INPUT the input path of the RMG-Java thermo library directory -LIBRARYNAME the libraryname for the RMG-Py format thermo library - -""" - -import argparse -import os -from rmgpy.data.thermo import ThermoLibrary -from rmgpy import settings - -if __name__ == '__main__': - - parser = argparse.ArgumentParser() - parser.add_argument('inputPath', metavar='INPUT', type=str, nargs=1, - help='the input path of the RMG-Java thermo library directory') - parser.add_argument('libraryName', metavar='OUTPUT', type=str, nargs=1, - help='the libraryName for the RMG-Py format thermo library') - - args = parser.parse_args() - inputPath = args.inputPath[0] - libraryName = args.libraryName[0] - - library = ThermoLibrary() - library.loadOld( - dictstr = os.path.join(inputPath, 'Dictionary.txt'), - treestr = '', - libstr = os.path.join(inputPath, 'Library.txt'), - numParameters = 12, - numLabels = 1, - pattern = False, - ) - library.name = libraryName - - # Save in Py format - library.save(os.path.join(settings['database.directory'], 'thermo', 'libraries', libraryName+'.py')) diff --git a/scripts/importOldDatabase.py b/scripts/importOldDatabase.py deleted file mode 100644 index 8364d09ef4..0000000000 --- a/scripts/importOldDatabase.py +++ /dev/null @@ -1,44 +0,0 @@ -#!/usr/bin/env python -# encoding: utf-8 - -""" -This script imports an RMG-Java database from the output directory and saves -it in the user-specified database directory. - -The outputDirectory can then be used to overwrite the existing RMG-database files. -""" - -import os -import argparse - -from rmgpy.data.rmg import RMGDatabase - -################################################################################ - - -if __name__ == '__main__': - - parser = argparse.ArgumentParser() - parser.add_argument('inputPath', metavar='INPUT', type=str, nargs=1, - help='the input path of the RMG-Java database directory') - parser.add_argument('outputPath', metavar='OUTPUT', type=str, nargs=1, - help='output path for the desired RMG-Py database directory') - - args = parser.parse_args() - inputPath = args.inputPath[0] - outputPath = args.outputPath[0] - - newPath = 'input' - - print 'Loading old RMG-Java database...' - database = RMGDatabase() - database.loadOld(inputPath) - - try: - os.makedirs(outputPath) - except: - pass - - print 'Saving the new RMG-Py database...' - database.save(outputPath) - print "Done!" diff --git a/scripts/kineticsGroups.py b/scripts/kineticsGroups.py deleted file mode 100644 index eab3923351..0000000000 --- a/scripts/kineticsGroups.py +++ /dev/null @@ -1,448 +0,0 @@ -#!/usr/bin/env python -# encoding: utf-8 - -""" -This script is used for working with the kinetics group additivity values in -RMG. There are several different types of operations this script can do, and -each of these has a number of required and optional command-line arguments. -Use the "-h" flag to get more information. -""" - -import argparse -import os.path -import time -import math -import numpy -import pylab -import scipy.stats -import matplotlib -matplotlib.rc('mathtext', fontset='stixsans', default='regular') -import re -import rmgpy -from rmgpy.quantity import constants -from rmgpy.kinetics import Arrhenius, ArrheniusEP, KineticsData -from rmgpy.data.base import getAllCombinations -from rmgpy.species import Species -from rmgpy import settings - -################################################################################ - -def loadDatabase(): - print 'Loading RMG database...' - from rmgpy.data.rmg import RMGDatabase - database = RMGDatabase() - database.load(settings['database.directory'], kineticsFamilies='all', kineticsDepositories='all') - return database - -def convertKineticsToPerSiteBasis(kinetics, degeneracy): - """ - Given high-pressure-limit `kinetics` which includes reaction-path - `degeneracy`, convert the kinetics to be on a per-site basis. - """ - if isinstance(kinetics, KineticsData): - kinetics.kdata.value_si /= degeneracy - elif isinstance(kinetics, Arrhenius): - kinetics.A.value_si /= degeneracy - elif isinstance(kinetics, ArrheniusEP): - kinetics.A.value_si /= degeneracy - else: - raise Exception('Unable to convert kinetics of type {0} to per-site basis.'.format(kinetics.__class__)) - return kinetics - -################################################################################ - -def createDataSet(label, family, database): - """ - Create a data set from the component of the kinetics `family` indicated by - the given `label`. The full RMG `database` must be loaded so that we can - get thermodynamics for some species. - """ - dataset = [] - - if label == 'rules': - for label, entries in family.rules.entries.items(): - for entry in entries: - # Skip ArrheniusEP entries with Evans-Polanyi values - if isinstance(entry.data, ArrheniusEP) and entry.data.alpha.value != 0: continue - # Also skip entries with rank of zero (since they are just made-up numbers) - if entry.rank == 0: continue - template = [family.groups.entries[node] for node in label.split(';')] - reaction = entry.item - dataset.append([reaction, template, entry]) - else: - label = '{0}/{1}'.format(family.label, label) - for depository in family.depositories: - if depository.label == label: - break - else: - raise ValueError('Invalid value "{0}" for label parameter.'.format(label)) - - for entry in depository.entries.values(): - reaction, template = database.kinetics.getForwardReactionForFamilyEntry(entry=entry, family=family.label, thermoDatabase=database.thermo) - dataset.append([reaction, template, entry]) - - return dataset - -################################################################################ - -class ArgumentError(Exception): - """ - An exception raised when the command-line arguments given to the script are - invalid. Pass a string describing why the arguments are invalid. - """ - pass - -################################################################################ - -def generate(args): - """ - Generate kinetics group additivity values for one (or more) reaction - families. The `args` parameter provides the results of parsing the - command-line arguments using argparse. - """ - # Make sure we have at least one family to generate values for - if len(args.family) == 0 and not args.all: - raise ArgumentError('No reaction families specified') - - # Make sure the method is valid - method = args.method - if method not in ['KineticsData', 'Arrhenius', 'Arrhenius2']: - raise ArgumentError('Invalid method "{0}" specified'.format(method)) - - # If training sets are not specified, 'training' and 'rules' are used - trainingSetLabels = args.training - if not trainingSetLabels: - trainingSetLabels = ['rules', 'training'] - - # Load the database - database = loadDatabase() - - # If --all flag was specified, use all reaction families - families = [] - if args.all: - families = database.kinetics.families.keys() - else: - families = args.family - - # Iterate over each family, generating and saving group values - for label in families: - family = database.kinetics.families[label] - family.addKineticsRulesFromTrainingSet(thermoDatabase=database.thermo) - - trainingSet = [] - for trainingSetLabel in trainingSetLabels: - for reaction, template, entry in createDataSet(trainingSetLabel, family, database): - kinetics = reaction.kinetics or entry.data - trainingSet.append((template, kinetics)) - - kunits = family.getRateCoefficientUnits() - - # Generate the group values (implemented on the KineticsGroups class) - changed = family.groups.generateGroupAdditivityValues(trainingSet, kunits, method=method) - - if changed: - # Save the new group values to disk - family.saveGroups(os.path.join(settings['database.directory'], 'kinetics', 'families', label, 'groups.py')) - -################################################################################ - -def evaluate(args): - """ - Evaluate kinetics group additivity values for one (or more) reaction - families. The `args` parameter provides the results of parsing the - command-line arguments using argparse. - """ - method = 'rate rules' if args.rules else 'group additivity' - plot = 'interactive' if args.interactive else 'normal' - - # If test sets are not specified, choose some - testSets = args.test - if not testSets: - testSets = ['NIST'] - - # Load the database - database = loadDatabase() - - # If --all flag was specified, use all reaction families - families = [] - if args.all: - families = database.kinetics.families.keys() - else: - families = args.family - - # Iterate over each family, generating and saving group values - for label in families: - family = database.kinetics.families[label] - family.addKineticsRulesFromTrainingSet(thermoDatabase=database.thermo) - changed = evaluateKineticsGroupValues( - database = database, - family = family, - method = method, - testSetLabels = testSets, - plot = plot, - exactOnly = args.exact, - estimateOnly = args.estimate, - ) - -def evaluateKineticsGroupValues(family, database, method, testSetLabels, plot, exactOnly=False, estimateOnly=False): - """ - Evaluate the kinetics group additivity values for the given reaction - `family` using the specified lists of depository components - `testSetLabels` as the test set. The already-loaded RMG database should be - given as the `database` parameter. - """ - kunits = family.getRateCoefficientUnits() - - assert not (exactOnly and estimateOnly) - - print 'Categorizing reactions in test sets for {0}'.format(family.label) - testSets0 = [] - for testSetLabel in testSetLabels: - testSet = createDataSet(testSetLabel, family, database) - testSets0.append((testSetLabel, testSet)) - - for testSetLabel, testSet in testSets0: - for index in range(len(testSet)): - reaction, template, entry = testSet[index] - for reactant in reaction.reactants: - if isinstance(reactant, Species) and not reactant.label and len(reactant.molecule) > 0: - reactant.label = reactant.molecule[0].toSMILES() - for product in reaction.products: - if isinstance(product, Species) and not product.label and len(product.molecule) > 0: - product.label = product.molecule[0].toSMILES() - - # For each entry in each test set, determine the kinetics as predicted by - # RMG-Py and as given by the entry in the test set - # Note that this is done on a per-site basis! - kineticsModels = []; kineticsData = [] - testSets = [] - for testSetLabel, testSet0 in testSets0: - testSet = [] - for index in range(len(testSet0)): - reaction, template, entry = testSet0[index] - krule = family.getKineticsForTemplate(template, degeneracy=1, method='rate rules') - kgroup = family.getKineticsForTemplate(template, degeneracy=1, method='group additivity') - kdata = convertKineticsToPerSiteBasis(reaction.kinetics, reaction.degeneracy) - if exactOnly and not re.search('Exact', krule.comment): - continue - elif estimateOnly and not re.search('Estimated', krule.comment): - continue - testSet.append((reaction, template, entry, krule, kgroup, kdata)) - testSets.append((testSetLabel, testSet)) - - # Generate parity plots at several temperatures - print 'Generating parity plots for {0}'.format(family.label) - - import matplotlib.pyplot as plt - from matplotlib.widgets import CheckButtons - - Tdata = [500,1000,1500,2000] - - if kunits == 'm^3/(mol*s)': - kunits = 'cm$^3$/mol*s'; kfactor = 1.0e6 - elif kunits == 's^-1': - kunits = 's$^{-1}$'; kfactor = 1.0 - - for T in Tdata: - - stdev_total = 0; ci_total = 0; count_total = 0 - - # Initialize plot - if plot == 'interactive': - fig = pylab.figure(figsize=(10,8)) - ax = plt.subplot(1, 1, 1) - else: - fig = pylab.figure(figsize=(6,5)) - ax = plt.subplot(1, 1, 1) - ax = plt.subplot(1, 1, 1) - lines = [] - legend = [] - - # Iterate through the test sets, plotting each - for testSetLabel, testSet in testSets: - - kmodel = []; kdata = [] - stdev = 0; ci = 0; count = 0 - - for reaction, template, entry, kineticsRule, kineticsGroup, kineticsData in testSet: - - if method == 'rate rules': - kineticsModel = kineticsRule - elif method == 'group additivity': - kineticsModel = kineticsGroup - - # Honor temperature ranges when plotting data - # Place a dummy value so that the points so that the - # interactivity is still correct - if not kineticsData.isTemperatureValid(T): - kmodel.append(0.0) - kdata.append(0.0) - continue - - # Evaluate k(T) for both model and data at this temperature - if isinstance(kineticsModel, ArrheniusEP): - km = kineticsModel.getRateCoefficient(T, 0) * kfactor - else: - km = kineticsModel.getRateCoefficient(T) * kfactor - kmodel.append(km) - if isinstance(kineticsData, ArrheniusEP): - kd = kineticsData.getRateCoefficient(T, 0) * kfactor - else: - kd = kineticsData.getRateCoefficient(T) * kfactor - kdata.append(kd) - - # Evaluate variance - stdev += (math.log10(km) - math.log10(kd))**2 - count += 1 - - stdev_total += stdev - count_total += count - stdev = math.sqrt(stdev / (count - 1)) - ci = scipy.stats.t.ppf(0.975, count - 1) * stdev - - assert len(kmodel) == len(testSet) - assert len(kdata) == len(testSet) - - print "Test set {0} contained {1} rates.".format(testSetLabel, count) - print 'Confidence interval at T = {0:g} K for test set "{1}" = 10^{2:g}'.format(T, testSetLabel, ci) - - # Add this test set to the plot - lines.append(ax.loglog(kdata, kmodel, 'o', picker=5)[0]) - legend.append(testSetLabel) - - stdev_total = math.sqrt(stdev_total / (count_total - 1)) - ci_total = scipy.stats.t.ppf(0.975, count_total - 1) * stdev_total - - print 'Total confidence interval at T = {0:g} K for all test sets = 10^{1:g}'.format(T, ci_total) - - # Finish plots - xlim = pylab.xlim() - ylim = pylab.ylim() - lim = (min(xlim[0], ylim[0])*0.1, max(xlim[1], ylim[1])*10) - ax.loglog(lim, lim, '-k') - ax.loglog(lim, [lim[0] * 10**ci_total, lim[1] * 10**ci_total], '--k') - ax.loglog(lim, [lim[0] / 10**ci_total, lim[1] / 10**ci_total], '--k') - pylab.xlabel('Actual rate coefficient ({0})'.format(kunits)) - pylab.ylabel('Predicted rate coefficient ({0})'.format(kunits)) - if len(testSets) > 1: - pylab.legend(legend, loc=4, numpoints=1) - pylab.title('%s, T = %g K' % (family.label, T)) - pylab.xlim(lim) - pylab.ylim(lim) - - plot_range = math.log10(lim[1] / lim[0]) - if plot_range > 25: - majorLocator = matplotlib.ticker.LogLocator(1e5) - minorLocator = matplotlib.ticker.LogLocator(1e5, subs=[1, 10, 100, 1000, 10000]) - elif plot_range > 20: - majorLocator = matplotlib.ticker.LogLocator(1e4) - minorLocator = matplotlib.ticker.LogLocator(1e4, subs=[1, 10, 100, 1000]) - elif plot_range > 15: - majorLocator = matplotlib.ticker.LogLocator(1e3) - minorLocator = matplotlib.ticker.LogLocator(1e3, subs=[1, 10, 100]) - elif plot_range > 10: - majorLocator = matplotlib.ticker.LogLocator(1e2) - minorLocator = matplotlib.ticker.LogLocator(1e2, subs=[1, 10]) - else: - majorLocator = matplotlib.ticker.LogLocator(1e1) - minorLocator = None - ax.xaxis.set_major_locator(majorLocator) - ax.yaxis.set_major_locator(majorLocator) - if minorLocator: - ax.xaxis.set_minor_locator(minorLocator) - ax.yaxis.set_minor_locator(minorLocator) - - def onpick(event): - index = lines.index(event.artist) - xdata = event.artist.get_xdata() - ydata = event.artist.get_ydata() - testSetLabel, testSet = testSets[index] - for ind in event.ind: - reaction, template, entry, krule, kgroup, kdata = testSet[ind] - kunits = 'm^3/(mol*s)' if len(reaction.reactants) == 2 else 's^-1' - print testSetLabel - print 'template = [{0}]'.format(', '.join([g.label for g in template])) - print 'entry = {0!r}'.format(entry) - print str(reaction) - print 'k_data = {0:9.2e} {1}'.format(xdata[ind], kunits) - print 'k_model = {0:9.2e} {1}'.format(ydata[ind], kunits) - print krule - if kgroup: print kgroup - print krule.comment - if kgroup: print kgroup.comment - print - - connection_id = fig.canvas.mpl_connect('pick_event', onpick) - - if plot == 'interactive': - rax = plt.axes([0.15, 0.65, 0.2, 0.2]) - check = CheckButtons(rax, legend, [True for label in legend]) - - def func(label): - for index in range(len(lines)): - if legend[index] == label: - lines[index].set_visible(not lines[index].get_visible()) - plt.draw() - check.on_clicked(func) - - fig.subplots_adjust(left=0.10, bottom=0.10, right=0.97, top=0.95, wspace=0.20, hspace=0.20) - - else: - fig.subplots_adjust(left=0.15, bottom=0.14, right=0.95, top=0.93, wspace=0.20, hspace=0.20) - filename = '{0}_{1:g}'.format(family.label, T) - if method == 'rate rules': - filename += '_rules' - elif method == 'group additivity': - filename += '_groups' - if exactOnly: - filename += '_exact' - elif estimateOnly: - filename += '_estimate' - pylab.savefig('{0}.pdf'.format(filename)) - pylab.savefig('{0}.png'.format(filename), dpi=200) - - pylab.show() - -################################################################################ - -def parseAndRunCommandLineArguments(): - - parser = argparse.ArgumentParser() - subparsers = parser.add_subparsers(dest='command', help='') - - # generate - generate and save kinetics group additivity values - generateParser = subparsers.add_parser('generate', help='generate and save kinetics group values for one or more families') - generateParser.add_argument('family', metavar='', type=str, nargs='*', help='the family to generate, or --all for all families') - generateParser.add_argument('-a', '--all', action='store_true', help='generate for all families') - generateParser.add_argument('-m', '--method', metavar='', type=str, nargs='?', default='Arrhenius', help='the method to use') - generateParser.add_argument('--training', metavar='', type=str, nargs='*', help='the training set(s) to use') - generateParser.set_defaults(run=generate) - - # evaluate - load and evaluate kinetics group additivity values - evaluateParser = subparsers.add_parser('evaluate', help='evaluate kinetics group values for one family') - evaluateParser.add_argument('family', metavar='', type=str, nargs=1, help='the family to evaluate') - evaluateParser.add_argument('-a', '--all', action='store_true', help='generate for all families') - evaluateParser.add_argument('--test', metavar='', type=str, nargs='*', help='the test set(s) to use') - evaluateParser.add_argument('-i', '--interactive', action='store_true', help='evaluate using interactive plots') - evaluateParser.add_argument('--rules', action='store_true', help='use rate rules instead of group additivity') - evaluateParser.add_argument('--exact', action='store_true', help='only plot exact matches') - evaluateParser.add_argument('--estimate', action='store_true', help='only plot estimated matches') - evaluateParser.set_defaults(run=evaluate) - - args = parser.parse_args() - - try: - args.run(args) - except ArgumentError, e: - for choice, subparser in subparsers.choices.iteritems(): - if args.command == choice: - subparser.print_help() - break - else: - parser.print_help() - print 'ArgumentError: {0}'.format(e) - -################################################################################ - -if __name__ == '__main__': - parseAndRunCommandLineArguments() diff --git a/scripts/kineticsTraining.py b/scripts/kineticsTraining.py deleted file mode 100644 index f9bccbd8de..0000000000 --- a/scripts/kineticsTraining.py +++ /dev/null @@ -1,487 +0,0 @@ -#!/usr/bin/env python -# encoding: utf-8 - -""" -This script can be used to generate the "best fit" high-pressure limit kinetics -to a set of input kinetics. - -To simply generate the best-fit kinetics, use the command :: - -$ python kineticsTraining.py generate - -where is the name of the reaction family and is the index of -the reaction to generate the recommended kinetics for. To also generate a plot -of the best-fit kinetics, use the command :: - -$ python kineticsTraining.py evaluate - -""" - -import os.path -import math -import numpy -import matplotlib -matplotlib.rc('mathtext', fontset='stixsans', default='regular') -import pylab - -from rmgpy.quantity import Quantity, constants -from rmgpy.thermo import * -from rmgpy.kinetics import * -from rmgpy.data.reference import * -from rmgpy.data.base import Entry -from rmgpy.data.thermo import ThermoDatabase -from rmgpy.data.kinetics import saveEntry -from rmgpy.molecule import Molecule -from rmgpy.species import Species -from rmgpy.reaction import Reaction - -################################################################################ - -class ArgumentError(Exception): - """ - An exception raised when the command-line arguments given to the script are - invalid. Pass a string describing why the arguments are invalid. - """ - pass - -################################################################################ - -forwardKinetics = [] -reverseKinetics = [] -forwardReaction = None -reverseReaction = None -forwardWeights = [] -reverseWeights = [] - -def loadSpecies(adjlist): - species = Species().fromAdjacencyList(adjlist) - species.molecule = species.molecule[0].generate_resonance_structures() - species.thermo = getThermoData(species) - species.molecule = [Molecule().fromAdjacencyList(adjlist)] - return species - -def reaction(index, label, reactant1, product1, forwardDegeneracy=1, reverseDegeneracy=1, reactant2=None, product2=None, product3=None): - global forwardReaction, reverseReaction - reactants = [loadSpecies(reactant1)] - if reactant2: - reactants.append(loadSpecies(reactant2)) - products = [loadSpecies(product1)] - if product2: - products.append(loadSpecies(product2)) - if product3: - reactants.append(loadSpecies(product3)) - forwardReaction = Reaction(reactants=reactants, products=products, degeneracy=forwardDegeneracy) - reverseReaction = Reaction(reactants=products, products=reactants, degeneracy=reverseDegeneracy) - -def entry(forward, kinetics, reference, referenceType, shortDesc, longDesc, weight=1.0): - global forwardKinetics, reverseKinetics - referenceTypes = {'theory': 'T', 'experiment': 'E', 'review': 'R'} - try: - comment = referenceTypes[referenceType] + '/' - except KeyError: - comment = '' - author = reference.authors[0].split() - if author[-1] == 'Jr.': - author = author[-2][:-1] - else: - author = author[-1] - comment += '{0!s}/{1!s}'.format(reference.year, author) - kinetics.comment = comment - if forward: - forwardKinetics.append(kinetics) - forwardWeights.append(weight) - else: - reverseKinetics.append(kinetics) - reverseWeights.append(weight) - -def loadKinetics(path): - """ - Load a set of kinetics data from the file located at `path` on disk. - """ - global forwardKinetics, reverseKinetics - global forwardReaction, reverseReaction - global forwardWeights, reverseWeights - - print 'Loading kinetics from {0}...'.format(os.path.relpath(path)) - - forwardKinetics = [] - reverseKinetics = [] - forwardReaction = None - reverseReaction = None - forwardWeights = [] - reverseWeights = [] - - # Set up global and local context - global_context = { - '__builtins__': None, - 'True': True, - 'False': False, - 'forwardKinetics': [], - 'reverseKinetics': [], - 'forwardWeights': [], - 'reverseWeights': [], - 'forwardReaction': None, - 'reverseReaction': None, - } - local_context = { - '__builtins__': None, - # Function prototypes - 'reaction': reaction, - 'entry': entry, - # Constants - 'R': constants.R, - 'kB': constants.kB, - # Kinetics types - 'KineticsData': KineticsData, - 'Arrhenius': Arrhenius, - 'ArrheniusEP': ArrheniusEP, - 'MultiKinetics': MultiKinetics, - 'PDepArrhenius': PDepArrhenius, - 'Chebyshev': Chebyshev, - 'ThirdBody': ThirdBody, - 'Lindemann': Lindemann, - 'Troe': Troe, - # Reference types - 'Reference': Reference, - 'Article': Article, - 'Book': Book, - 'Thesis': Thesis, - } - - # Process the file - f = open(path, 'r') - exec f in global_context, local_context - f.close() - - # For each kinetics entry in reverseKinetics, fit the reverse kinetics and - # append it to forwardKinetics - for kinetics0, weight in zip(reverseKinetics, reverseWeights): - reverseReaction.kinetics = kinetics0 - try: - kinetics = reverseReaction.generateReverseRateCoefficient() - except IndexError: - continue - kinetics.Tmin = kinetics0.Tmin - kinetics.Tmax = kinetics0.Tmax - kinetics.comment = kinetics0.comment + '*' - forwardKinetics.append(kinetics) - forwardWeights.append(weight) - - return forwardKinetics, reverseKinetics, forwardWeights, reverseWeights - -################################################################################ - -thermoDatabase = None - -def loadThermoDatabase(path): - """ - Load the RMG thermodynamics database from `path`. - """ - global thermoDatabase - print 'Loading thermodynamics database...' - thermoDatabase = ThermoDatabase() - thermoDatabase.load(path) - -def getThermoData(species): - global thermoDatabase - - # Ensure molecules are using explicit hydrogens - implicitH = [mol.implicitHydrogens for mol in species.molecule] - for molecule in species.molecule: - molecule.makeHydrogensExplicit() - - thermo = [] - for molecule in species.molecule: - molecule.clearLabeledAtoms() - molecule.updateAtomTypes() - thermo.append(thermoDatabase.getThermoData(species)) - - H298 = numpy.array([t.getEnthalpy(298.) for t in thermo]) - indices = H298.argsort() - - # If multiple resonance isomers are present, use the thermo data of - # the most stable isomer (i.e. one with lowest enthalpy of formation) - # as the thermo data of the species - thermo0 = thermo[indices[0]] - - # Sort the structures in order of decreasing stability - species.molecule = [species.molecule[ind] for ind in indices] - implicitH = [implicitH[ind] for ind in indices] - - # Convert to Wilhoit - if isinstance(thermo0, Wilhoit): - wilhoit = thermo0 - else: - linear = species.molecule[0].isLinear() - nRotors = species.molecule[0].countInternalRotors() - nFreq = 3 * len(species.molecule[0].atoms) - (5 if linear else 6) - nRotors - wilhoit = convertThermoModel(thermo0, Wilhoit, linear=linear, nFreq=nFreq, nRotors=nRotors) - - # Restore implicit hydrogens if necessary - for implicit, molecule in zip(implicitH, species.molecule): - if implicit: molecule.makeHydrogensImplicit() - - return wilhoit - -################################################################################ - -def fitArrhenius(kineticsList, weights, Tlist, Tmin, Tmax): - """ - Fit a modified Arrhenius expression to the set of loaded kinetics - `kineticsList` using the given array of temperatures `Tlist` in K. The - `Tmin` and `Tmax` parameters specify the limits in K at which the fitted - kinetics should be said to be valid. Returns the best-fit kinetics and the - confidence interval (on a log scale). - """ - import scipy.stats - - Tdata = []; kdata = []; wdata = [] - - print 'Fitting modified Arrhenius kinetics...' - - for kinetics, weight in zip(kineticsList, weights): - for n in range(len(Tlist)): - if kinetics.isTemperatureValid(Tlist[n]): - Tdata.append(Tlist[n]) - kdata.append(kinetics.getRateCoefficient(Tlist[n])) - wdata.append(weight) - - Tdata = numpy.array(Tdata, numpy.float) - kdata = numpy.array(kdata, numpy.float) - wdata = numpy.array(wdata, numpy.float) - arrhenius = Arrhenius().fitToData(Tdata, kdata, kunits='m^3/(mol*s)', T0=300, weights=wdata) - - # Compute RMS error and confidence interval - count = len(kdata) - rmsError = 0 - for T, k in zip(Tdata, kdata): - rmsError += (math.log10(k) - math.log10(arrhenius.getRateCoefficient(T)))**2 - rmsError = math.sqrt(rmsError / (count - 3)) - ci = scipy.stats.t.ppf(0.975, count - 3) * rmsError - - # Adjust units of best-fit Arrhenius expression - arrhenius.A.units = 'cm^3/(mol*s)' if len(forwardReaction.reactants) == 2 else 's^-1' - arrhenius.Ea.units = 'kJ/mol' - - # Set Tmin and Tmax of best-fit Arrhenius expression - arrhenius.Tmin = Quantity(Tmin,"K") - arrhenius.Tmax = Quantity(Tmax,"K") - - return arrhenius, ci - -################################################################################ - -def plotKinetics(kineticsList, Tlist, filename=None, bestKinetics=None): - """ - Plot the set of loaded kinetics `kineticsList` at the array of temperatures - `Tlist` in K. Different symbols denote the various reference types, while - different linespecs denote individual kinetics within each reference type. - If given, the `bestKinetics` will also be plotted, using a thicker line to - make it stand out. - """ - - sm = pylab.cm.ScalarMappable( - norm = matplotlib.colors.Normalize(vmin=0, vmax=len(kineticsList)-1), - cmap = pylab.get_cmap('jet'), - ) - - if len(forwardReaction.reactants) == 2: - kfactor = 1e6; kunits = '$cm^3/mol*s$' - else: - kfactor = 1; kunits = '$s^-1$' - - fig = pylab.figure(figsize=(8,6)) - legend = []; lines = [] - for index, kinetics in enumerate(kineticsList): - klist = numpy.zeros_like(Tlist) - for n in range(len(Tlist)): - if kinetics.isTemperatureValid(Tlist[n]): - klist[n] = kinetics.getRateCoefficient(Tlist[n]) - try: - if kinetics.comment[0] == 'R': - linespec = '-' - elif kinetics.comment[0] == 'E': - linespec = ':' - elif kinetics.comment[0] == 'T': - linespec = '--' - else: - linespec = '-.' - except IndexError: - continue - color = sm.to_rgba(index) - if numpy.any(klist): - lines.append(pylab.semilogy(1000. / Tlist, klist * kfactor, linespec, color=color, picker=5)[0]) - legend.append(kinetics.comment) - - if bestKinetics: - klist = numpy.zeros_like(Tlist) - for n in range(len(Tlist)): - if bestKinetics.isTemperatureValid(Tlist[n]): - klist[n] = bestKinetics.getRateCoefficient(Tlist[n]) - if numpy.any(klist): - pylab.semilogy(1000. / Tlist, klist * kfactor, '-k', linewidth=3) - legend.append(bestKinetics.comment) - - pylab.xlabel('1000 / (Temperature (K))') - pylab.ylabel('Rate coefficient ({0})'.format(kunits)) - pylab.xlim(0.0,4.0) - pylab.legend(legend, loc=1) - - pylab.title('{0} $\\rightarrow$ {1}'.format( - ' + '.join([spec.label for spec in forwardReaction.reactants]), - ' + '.join([spec.label for spec in forwardReaction.products]), - )) - - if filename: - pylab.savefig(filename) - - def onpick(event): - index = lines.index(event.artist) - kinetics = kineticsList[index] - print kinetics.comment - print kinetics - - connection_id = fig.canvas.mpl_connect('pick_event', onpick) - - pylab.show() - -################################################################################ - -def generateEntry(reaction, kinetics, index, ci): - """ - Return a string containing the reaction and its recommended kinetics in a - format suitable for putting in the training set. - """ - import StringIO - - longDesc = "" - - entry = Entry( - index = index, - item = forwardReaction, - data = kinetics, - reference = None, - referenceType = "review", - shortDesc = "Recommended value based on evaluation of {0:d} kinetics entries. log(CI) = {1:.3f}".format(len(forwardKinetics), ci), - longDesc = longDesc, - ) - - f = StringIO.StringIO() - saveEntry(f, entry) - string = f.getvalue() - f.close() - - return string - -################################################################################ - -Tlist_fit = 1000.0/numpy.arange(0.5, 3.34, 0.01, numpy.float) -Tlist_plot = 1000.0/numpy.arange(0.5, 3.34, 0.1, numpy.float) -Tmin = 300 -Tmax = 2000 - -def getPath(family, index): - from rmgpy import settings - return os.path.join(settings['database.directory'], 'kinetics', 'families', family, 'training', '{0}.py'.format(index)) - -def generate(args): - """ - Generate the recommended kinetics for a reaction. - """ - family = str(args.family[0]) - - for index in args.index: - index = int(index) - - path = getPath(args.family[0], str(index)) - forwardKinetics, reverseKinetics, forwardWeights, reverseWeights = loadKinetics(path) - kinetics, ci = fitArrhenius(forwardKinetics, forwardWeights, Tlist_fit, Tmin, Tmax) - kinetics.comment = 'Best fit' - - print - print generateEntry(forwardReaction, kinetics, index, ci) - -def evaluate(args): - """ - Evaluate the collected kinetics for a given reaction. - """ - family = str(args.family[0]) - - for index in args.index: - index = int(index) - - path = getPath(args.family[0], str(index)) - forwardKinetics, reverseKinetics, forwardWeights, reverseWeights = loadKinetics(path) - - done = False - while not done: - kinetics, ci = fitArrhenius(forwardKinetics, forwardWeights, Tlist_fit, Tmin, Tmax) - kinetics.comment = 'Best fit' - - Tlist = Tlist_fit - klist = numpy.zeros_like(Tlist) - for n in range(len(Tlist)): - if kinetics.isTemperatureValid(Tlist[n]): - klist[n] = kinetics.getRateCoefficient(Tlist[n]) - - done = True; outlier = None; outlierRMS = 0 - for kinetics0 in forwardKinetics: - klist0 = numpy.zeros_like(Tlist) - rms = 0; count = 0 - for n in range(len(Tlist)): - if kinetics0.isTemperatureValid(Tlist[n]): - klist0[n] = kinetics0.getRateCoefficient(Tlist[n]) - error = math.log10(klist0[n]) - math.log10(klist[n]) - rms += error * error - count += 1 - if count == 0: continue - rms = math.sqrt(rms / count) - if rms > 3.0 and rms > outlierRMS: - outlier = kinetics0 - outlierRMS = rms - if outlier is not None: - forwardKinetics.remove(outlier) - print 'Identified kinetics "{0}" as outlier (RMS = {1:g}). Removing it from Arrhenius fit.'.format(outlier.comment, outlierRMS) - done = False - - plotKinetics(forwardKinetics, Tlist_plot, filename='{0}.pdf'.format(path[:-3]), bestKinetics=kinetics) - - print - print generateEntry(forwardReaction, kinetics, index, ci) - -################################################################################ - -if __name__ == '__main__': - - import argparse - from rmgpy import settings - - parser = argparse.ArgumentParser() - subparsers = parser.add_subparsers(dest='command', help='') - - # generate - generate the recommended kinetics for a given reaction - generateParser = subparsers.add_parser('generate', help='generate the recommended kinetics for a given reaction') - generateParser.add_argument('family', type=str, nargs=1, help='the reaction family containing the reaction') - generateParser.add_argument('index', type=str, nargs='+', help='the index of the reaction') - generateParser.set_defaults(run=generate) - - # evaluate - evaluate the collected kinetics for a given reaction - evaluateParser = subparsers.add_parser('evaluate', help='evaluate the collected kinetics for a given reaction') - evaluateParser.add_argument('family', type=str, nargs=1, help='the reaction family containing the reaction') - evaluateParser.add_argument('index', metavar='', type=str, nargs='+', help='the index of the reaction') - evaluateParser.set_defaults(run=evaluate) - - args = parser.parse_args() - - loadThermoDatabase(os.path.join(settings['database.directory'], 'thermo')) - - try: - args.run(args) - except ArgumentError, e: - for choice, subparser in subparsers.choices.iteritems(): - if args.command == choice: - subparser.print_help() - break - else: - parser.print_help() - print 'ArgumentError: {0}'.format(e) diff --git a/scripts/showNewTrainingRxnsImprovements.ipynb b/scripts/showNewTrainingRxnsImprovements.ipynb index 20875c2714..fc289ce508 100644 --- a/scripts/showNewTrainingRxnsImprovements.ipynb +++ b/scripts/showNewTrainingRxnsImprovements.ipynb @@ -11,17 +11,15 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "from rmgpy.data.rmg import RMGDatabase\n", - "from rmgpy.chemkin import saveChemkinFile, saveSpeciesDictionary\n", - "from rmgpy.rmg.model import Species, getFamilyLibraryObject, CoreEdgeReactionModel\n", + "from rmgpy.chemkin import save_chemkin_file, save_species_dictionary\n", + "from rmgpy.rmg.model import Species, get_family_library_object, CoreEdgeReactionModel\n", "from rmgpy import settings\n", "from IPython.display import display\n", - "from rmgpy.cantherm.output import prettify\n", + "from arkane.output import prettify\n", "from rmgpy.kinetics.kineticsdata import KineticsData\n", "from rmgpy.data.kinetics.family import TemplateReaction\n", "from rmgpy.data.kinetics.depository import DepositoryReaction\n", @@ -31,14 +29,13 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "database = RMGDatabase()\n", "libraries = ['vinylCPD_H']\n", - "database.load(settings['database.directory'], kineticsFamilies='all', reactionLibraries = libraries, kineticsDepositories='all')" + "database.load(settings['database.directory'], kinetics_families='all', reaction_libraries = libraries, \n", + " kinetics_depositories='all')" ] }, { @@ -51,15 +48,13 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ - "reactionDict = {}\n", - "for libraryName in libraries:\n", - " kineticLibrary = database.kinetics.libraries[libraryName]\n", - " for index, entry in kineticLibrary.entries.iteritems():\n", + "reaction_dict = {}\n", + "for library_name in libraries:\n", + " kinetic_library = database.kinetics.libraries[library_name]\n", + " for index, entry in kinetic_library.entries.items():\n", " lib_rxn = entry.item\n", " display(lib_rxn)\n", " # Let's make RMG generate this reaction from the families!\n", @@ -75,7 +70,7 @@ " for mutation_i in range(rxt_mol_mutation_num):\n", " rxts_mol = [spc.molecule[mutation_i%(len(spc.molecule))] for spc in lib_rxn.reactants]\n", " pdts_mol = [spc.molecule[0] for spc in lib_rxn.products]\n", - " fam_rxn_list.extend(database.kinetics.generateReactionsFromFamilies(\n", + " fam_rxn_list.extend(database.kinetics.generate_reactions_from_families(\n", " reactants=rxts_mol, products=pdts_mol))\n", "\n", " if len(fam_rxn_list) == 1:\n", @@ -84,7 +79,7 @@ " fam_reactants = [r for r in fam_rxn.reactants]\n", " for lib_reactant in lib_reactants:\n", " for fam_reactant in fam_reactants:\n", - " if lib_reactant.isIsomorphic(fam_reactant):\n", + " if lib_reactant.is_isomorphic(fam_reactant):\n", " fam_reactants.remove(fam_reactant)\n", " break\n", "\n", @@ -92,20 +87,20 @@ " fam_products = [r for r in fam_rxn.products]\n", " for lib_product in lib_products:\n", " for fam_product in fam_products:\n", - " if lib_product.isIsomorphic(fam_product):\n", + " if lib_product.is_isomorphic(fam_product):\n", " fam_products.remove(fam_product)\n", " break\n", "\n", " forward = not (len(fam_reactants) != 0 or len(fam_products) != 0)\n", " # find the labeled atoms using family and reactants & products from fam_rxn \n", " family_database = database.kinetics.families[fam_rxn.family]\n", - " family_database.addAtomLabelsForReaction(fam_rxn)\n", + " family_database.add_atom_labels_for_reaction(fam_rxn)\n", " fam_rxn.index = lib_rxn.index\n", - " reactionDict[lib_rxn]=fam_rxn\n", + " reaction_dict[lib_rxn]=fam_rxn\n", " for spec in fam_rxn.reactants + fam_rxn.products:\n", - " print spec.molecule[0].toSMILES()\n", + " print(spec.molecule[0].to_smiles())\n", " else:\n", - " print 'There was an issue finding a single reaction family for this reaction. Please investigate.'" + " print('There was an issue finding a single reaction family for this reaction. Please investigate.')" ] }, { @@ -125,9 +120,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "# Pick a temperature to evaluate the kinetics\n", @@ -137,106 +130,118 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "cem = CoreEdgeReactionModel()\n", - "cem.kineticsEstimator = 'rate rules'\n", - "cem.verboseComments = True" + "cem.kinetics_estimator = 'rate rules'\n", + "cem.verbose_comments = True" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ - "for libraryName in libraries:\n", - " kineticLibrary = database.kinetics.libraries[libraryName]\n", - " for index, entry in kineticLibrary.entries.iteritems():\n", + "for library_name in libraries:\n", + " kinetic_library = database.kinetics.libraries[library_name]\n", + " for index, entry in kinetic_library.entries.items():\n", " lib_rxn = entry.item\n", " try:\n", - " fam_rxn = reactionDict[lib_rxn]\n", + " fam_rxn = reaction_dict[lib_rxn]\n", " except:\n", - " print 'There was an issue finding a single reaction family for this reaction in step 1. Skipping.'\n", + " print('There was an issue finding a single reaction family for this reaction in step 1. Skipping.')\n", " continue\n", " for idx, spec in enumerate(fam_rxn.reactants):\n", " spec = Species(label=spec.label, molecule=spec.molecule)\n", - " spec.getThermoData()\n", + " spec.get_thermo_data()\n", " fam_rxn.reactants[idx] = spec\n", " for idx, spec in enumerate(fam_rxn.products):\n", " spec = Species(label=spec.label, molecule=spec.molecule)\n", - " spec.getThermoData()\n", + " spec.get_thermo_data()\n", " fam_rxn.products[idx] = spec\n", "\n", - " family = getFamilyLibraryObject(fam_rxn.family)\n", + " family = get_family_library_object(fam_rxn.family)\n", "\n", " # Set the reaction kinetics\n", - " kinetics, source, entryFamily, isForward = cem.generateKinetics(fam_rxn)\n", + " kinetics, source, entry_family, is_forward = cem.generate_kinetics(fam_rxn)\n", " fam_rxn.kinetics = kinetics\n", " # Flip the reaction direction if the kinetics are defined in the reverse direction\n", - " if not isForward:\n", + " if not is_forward:\n", " fam_rxn.reactants, fam_rxn.products = fam_rxn.products, fam_rxn.reactants\n", " fam_rxn.pairs = [(p,r) for r,p in fam_rxn.pairs]\n", - " if family.ownReverse and hasattr(fam_rxn,'reverse'):\n", - " if not isForward:\n", + " if family.own_reverse and hasattr(fam_rxn,'reverse'):\n", + " if not is_forward:\n", " fam_rxn.template = fam_rxn.reverse.template\n", " # We're done with the \"reverse\" attribute, so delete it to save a bit of memory\n", " delattr(fam_rxn,'reverse')\n", "\n", " # convert KineticsData to Arrhenius forms\n", " if isinstance(fam_rxn.kinetics, KineticsData):\n", - " fam_rxn.kinetics = fam_rxn.kinetics.toArrhenius()\n", + " fam_rxn.kinetics = fam_rxn.kinetics.to_arrhenius()\n", " # correct barrier heights of estimated kinetics\n", " if isinstance(fam_rxn,TemplateReaction) or isinstance(fam_rxn,DepositoryReaction): # i.e. not LibraryReaction\n", - " fam_rxn.fixBarrierHeight() # also converts ArrheniusEP to Arrhenius.\n", + " fam_rxn.fix_barrier_height() # also converts ArrheniusEP to Arrhenius.\n", "\n", - " if cem.pressureDependence and fam_rxn.isUnimolecular():\n", + " if cem.pressure_dependence and fam_rxn.is_unimolecular():\n", " # If this is going to be run through pressure dependence code,\n", " # we need to make sure the barrier is positive.\n", - " fam_rxn.fixBarrierHeight(forcePositive=True)\n", - " print '==============='\n", - " print index\n", + " fam_rxn.fix_barrier_height(force_positive=True)\n", + " print('===============')\n", + " print(index)\n", " display(lib_rxn)\n", - " print ''\n", - " print 'Library Kinetics'\n", - " print prettify(repr(entry.data))\n", + " print('')\n", + " print('Library Kinetics')\n", + " print(prettify(repr(entry.data)))\n", " \n", - " print ''\n", - " print 'Reaction Family Kinetics After Training'\n", - " print 'Family: {}'.format(fam_rxn.family)\n", - " print prettify(repr(fam_rxn.kinetics))\n", + " print('')\n", + " print('Reaction Family Kinetics After Training')\n", + " print('Family: {}'.format(fam_rxn.family))\n", + " print(prettify(repr(fam_rxn.kinetics)))\n", " \n", - " print ''\n", - " print 'Kinetics evaluated at {} K'.format(T)\n", - " print 'Library Kinetics: {:.2E}'.format(entry.data.getRateCoefficient(T))\n", - " print 'Reaction Family Kinetics After Training: {:.2E}'.format(fam_rxn.kinetics.getRateCoefficient(T))" + " print('')\n", + " print('Kinetics evaluated at {} K'.format(T))\n", + " print('Library Kinetics: {:.2E}'.format(entry.data.get_rate_coefficient(T)))\n", + " print('Reaction Family Kinetics After Training: {:.2E}'.format(fam_rxn.kinetics.get_rate_coefficient(T)))" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python 2", + "display_name": "Python 3", "language": "python", - "name": "python2" + "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", - "version": 2 + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.11" + "pygments_lexer": "ipython3", + "version": "3.7.4" + }, + "pycharm": { + "stem_cell": { + "cell_type": "raw", + "source": [], + "metadata": { + "collapsed": false + } + } } }, "nbformat": 4, - "nbformat_minor": 0 -} + "nbformat_minor": 1 +} \ No newline at end of file