diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index d28ee47ce0..57c0ac353d 100644 --- a/rmgpy/chemkin.pyx +++ b/rmgpy/chemkin.pyx @@ -1209,18 +1209,16 @@ def read_species_block(f, species_dict, species_aliases, species_list): if token_upper == 'END': break - site_token = token.split('/')[0] - if site_token.upper() == 'SDEN': + species_name = token.split('/')[0] # CHO*/2/ indicates an adsorbate CHO* taking 2 surface sites + if species_name.upper() == 'SDEN': # SDEN/4.1e-9/ indicates surface site density continue # TODO actually read in the site density processed_tokens.append(token) - if token in species_dict: - species = species_dict[token] - elif site_token in species_dict: - species = species_dict[site_token] + if species_name in species_dict: + species = species_dict[species_name] else: - species = Species(label=token) - species_dict[token] = species + species = Species(label=species_name) + species_dict[species_name] = species species_list.append(species) diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx index f24defb8cc..017b26ff70 100644 --- a/rmgpy/kinetics/arrhenius.pyx +++ b/rmgpy/kinetics/arrhenius.pyx @@ -707,12 +707,52 @@ cdef class ArrheniusBM(KineticsModel): """ self._A.value_si *= factor + def to_cantera_kinetics(self): + """ + Converts the RMG ArrheniusBM object to a cantera BlowersMaselRate. + + BlowersMaselRate(A, b, Ea, W) where A is in units of m^3/kmol/s, + b is dimensionless, and Ea and W are in J/kmol + """ + import cantera as ct + + rate_units_conversion = {'1/s': 1, + 's^-1': 1, + 'm^3/(mol*s)': 1000, + 'm^6/(mol^2*s)': 1000000, + 'cm^3/(mol*s)': 1000, + 'cm^6/(mol^2*s)': 1000000, + 'm^3/(molecule*s)': 1000, + 'm^6/(molecule^2*s)': 1000000, + 'cm^3/(molecule*s)': 1000, + 'cm^6/(molecule^2*s)': 1000000, + } + + if self._T0.value_si != 1: + A = self._A.value_si / (self._T0.value_si) ** self._n.value_si + else: + A = self._A.value_si + + try: + A *= rate_units_conversion[self._A.units] # convert from /mol to /kmol + except KeyError: + raise ValueError(f'ArrheniusBM A-factor units {self._A.units} not found among accepted ' + 'units for converting to Cantera BlowersMaselRate object.') + + b = self._n.value_si + Ea = self._E0.value_si * 1000 # convert from J/mol to J/kmol + w = self._w0.value_si * 1000 # convert from J/mol to J/kmol + + return ct.BlowersMaselRate(A, b, Ea, w) + def set_cantera_kinetics(self, ct_reaction, species_list): """ - Sets a cantera Reaction() object with the modified Arrhenius object - converted to an Arrhenius form. + Accepts a cantera Reaction object and sets its rate to a Cantera BlowersMaselRate object. """ - raise NotImplementedError('set_cantera_kinetics() is not implemented for ArrheniusBM class kinetics.') + import cantera as ct + if not isinstance(ct_reaction.rate, ct.BlowersMaselRate): + raise TypeError("ct_reaction must have a cantera BlowersMaselRate as the rate attribute") + ct_reaction.rate = self.to_cantera_kinetics() ################################################################################ diff --git a/rmgpy/kinetics/surface.pyx b/rmgpy/kinetics/surface.pyx index 120cfc006d..708b9ac546 100644 --- a/rmgpy/kinetics/surface.pyx +++ b/rmgpy/kinetics/surface.pyx @@ -50,7 +50,7 @@ cdef class StickingCoefficient(KineticsModel): ======================= ============================================================= Attribute Description ======================= ============================================================= - `A` The preexponential factor + `A` The preexponential factor. Unitless (sticking probability) `T0` The reference temperature `n` The temperature exponent `Ea` The activation energy @@ -285,6 +285,34 @@ cdef class StickingCoefficient(KineticsModel): return string + def to_cantera_kinetics(self): + """ + Converts the RMG StickingCoefficient object to a cantera StickingArrheniusRate + + StickingArrheniusRate(A,b,E) where A is dimensionless, b is dimensionless, and E is in J/kmol + """ + import cantera as ct + import rmgpy.quantity + if type(self._A) != rmgpy.quantity.ScalarQuantity: + raise TypeError("A factor must be a dimensionless ScalarQuantity") + A = self._A.value_si + b = self._n.value_si + E = self._Ea.value_si * 1000 # convert from J/mol to J/kmol + + return ct.StickingArrheniusRate(A, b, E) + + + def set_cantera_kinetics(self, ct_reaction, species_list): + """ + Passes in a cantera Reaction() object and sets its + rate to a Cantera ArrheniusRate object. + """ + import cantera as ct + assert isinstance(ct_reaction.rate, ct.StickingArrheniusRate), "Must have a Cantera StickingArrheniusRate attribute" + + # Set the rate parameter to a cantera Arrhenius object + ct_reaction.rate = self.to_cantera_kinetics() + ################################################################################ cdef class StickingCoefficientBEP(KineticsModel): """ @@ -352,7 +380,7 @@ cdef class StickingCoefficientBEP(KineticsModel): self.Pmin, self.Pmax, self.coverage_dependence, self.comment)) property A: - """The preexponential factor.""" + """The preexponential factor. Dimensionless (sticking probability)""" def __get__(self): return self._A def __set__(self, value): @@ -395,7 +423,8 @@ cdef class StickingCoefficientBEP(KineticsModel): cpdef double get_sticking_coefficient(self, double T, double dHrxn=0.0) except -1: """ Return the sticking coefficient (dimensionless) at - temperature `T` in K and enthalpy of reaction `dHrxn` in J/mol. + temperature `T` in K and enthalpy of reaction `dHrxn` in J/mol. + Never exceeds 1.0. """ cdef double A, n, Ea, stickingCoefficient Ea = self.get_activation_energy(dHrxn) @@ -524,7 +553,7 @@ cdef class SurfaceArrhenius(Arrhenius): property A: """The preexponential factor. - This is the only thing different from a normal Arrhenius class.""" + This (and the coverage dependence) is the only thing different from a normal Arrhenius class.""" def __get__(self): return self._A def __set__(self, value): @@ -570,6 +599,56 @@ cdef class SurfaceArrhenius(Arrhenius): return (SurfaceArrhenius, (self.A, self.n, self.Ea, self.T0, self.Tmin, self.Tmax, self.Pmin, self.Pmax, self.coverage_dependence, self.uncertainty, self.comment)) + def to_cantera_kinetics(self): + """ + Converts the RMG SurfaceArrhenius object to a cantera InterfaceArrheniusRate + + InterfaceArrheniusRate(A,b,E) where A is in units like m^2/kmol/s (depending on dimensionality) + b is dimensionless, and E is in J/kmol + """ + import cantera as ct + + rate_units_conversion = {'1/s': 1, + 's^-1': 1, + 'm^2/(mol*s)': 1000, + 'm^4/(mol^2*s)': 1000000, + 'cm^2/(mol*s)': 1000, + 'cm^4/(mol^2*s)': 1000000, + 'm^2/(molecule*s)': 1000, + 'm^4/(molecule^2*s)': 1000000, + 'cm^2/(molecule*s)': 1000, + 'cm^4/(molecule^2*s)': 1000000, + 'cm^5/(mol^2*s)': 1000000, + 'm^5/(mol^2*s)': 1000000, + } + + if self._T0.value_si != 1: + A = self._A.value_si / (self._T0.value_si) ** self._n.value_si + else: + A = self._A.value_si + + try: + A *= rate_units_conversion[self._A.units] # convert from /mol to /kmol + except KeyError: + raise ValueError('Arrhenius A-factor units {0} not found among accepted units for converting to ' + 'Cantera Arrhenius object.'.format(self._A.units)) + + b = self._n.value_si + E = self._Ea.value_si * 1000 # convert from J/mol to J/kmol + return ct.InterfaceArrheniusRate(A, b, E) + + def set_cantera_kinetics(self, ct_reaction, species_list): + """ + Takes in a cantera Reaction object and sets its + rate to a cantera InterfaceArrheniusRate object. + """ + import cantera as ct + if not isinstance(ct_reaction.rate, ct.InterfaceArrheniusRate): + raise TypeError("ct_reaction.rate must be an InterfaceArrheniusRate") + + # Set the rate parameter to a cantera Arrhenius object + ct_reaction.rate = self.to_cantera_kinetics() + ################################################################################ cdef class SurfaceArrheniusBEP(ArrheniusEP): diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index f10e1dcefa..4407050988 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -245,6 +245,7 @@ def to_chemkin(self, species_list=None, kinetics=True): else: return rmgpy.chemkin.write_reaction_string(self) + def to_cantera(self, species_list=None, use_chemkin_identifier=False): """ Converts the RMG Reaction object to a Cantera Reaction object @@ -288,7 +289,10 @@ def to_cantera(self, species_list=None, use_chemkin_identifier=False): if self.kinetics: if isinstance(self.kinetics, Arrhenius): # Create an Elementary Reaction - ct_reaction = ct.Reaction(reactants=ct_reactants, products=ct_products, rate=ct.ArrheniusRate()) + if isinstance(self.kinetics, SurfaceArrhenius): # SurfaceArrhenius inherits from Arrhenius + ct_reaction = ct.Reaction(reactants=ct_reactants, products=ct_products, rate=ct.InterfaceArrheniusRate()) + else: + ct_reaction = ct.Reaction(reactants=ct_reactants, products=ct_products, rate=ct.ArrheniusRate()) elif isinstance(self.kinetics, MultiArrhenius): # Return a list of elementary reactions which are duplicates ct_reaction = [ct.Reaction(reactants=ct_reactants, products=ct_products, rate=ct.ArrheniusRate()) @@ -339,6 +343,14 @@ def to_cantera(self, species_list=None, use_chemkin_identifier=False): reactants=ct_reactants, products=ct_products, rate=rate ) + elif isinstance(self.kinetics, SurfaceArrhenius): + rate = self.kinetics.to_cantera_kinetics() + ct_reaction = ct.InterfaceReaction(equation=str(self), rate=rate) + + elif isinstance(self.kinetics, StickingCoefficient): + rate = self.kinetics.to_cantera_kinetics() + ct_reaction = ct.Reaction(equation=str(self), rate=rate) + elif isinstance(self.kinetics, Lindemann): high_rate = self.kinetics.arrheniusHigh.to_cantera_kinetics(arrhenius_class=True) low_rate = self.kinetics.arrheniusLow.to_cantera_kinetics(arrhenius_class=True) @@ -356,6 +368,9 @@ def to_cantera(self, species_list=None, use_chemkin_identifier=False): reactants=ct_reactants, products=ct_products, rate=rate ) + elif isinstance(self.kinetics, StickingCoefficient): + ct_reaction = ct.Reaction(reactants=ct_reactants, products=ct_products, rate=ct.StickingArrheniusRate()) + else: raise NotImplementedError('Unable to set cantera kinetics for {0}'.format(self.kinetics)) @@ -1174,7 +1189,7 @@ def is_balanced(self): from rmgpy.molecule.element import element_list from rmgpy.molecule.fragment import CuttingLabel, Fragment - cython.declare(reactant_elements=dict, product_elements=dict, molecule=Graph, atom=Vertex, element=Element) + cython.declare(reactant_elements=dict, product_elements=dict, molecule=Molecule, atom=Atom, element=Element) reactant_elements = {} product_elements = {} diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index 8df771980d..e1a799da4e 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -79,7 +79,8 @@ from rmgpy.thermo.thermoengine import submit from rmgpy.tools.plot import plot_sensitivity from rmgpy.tools.uncertainty import Uncertainty, process_local_results -from rmgpy.yml import RMSWriter +from rmgpy.yaml_rms import RMSWriter +from rmgpy.yaml_cantera import CanteraWriter from rmgpy.rmg.reactors import Reactor ################################################################################ @@ -724,6 +725,7 @@ def register_listeners(self): self.attach(ChemkinWriter(self.output_directory)) self.attach(RMSWriter(self.output_directory)) + self.attach(CanteraWriter(self.output_directory)) if self.generate_output_html: self.attach(OutputHTMLWriter(self.output_directory)) @@ -1755,7 +1757,7 @@ def generate_cantera_files(self, chemkin_file, **kwargs): """ transport_file = os.path.join(os.path.dirname(chemkin_file), "tran.dat") file_name = os.path.splitext(os.path.basename(chemkin_file))[0] + ".yaml" - out_name = os.path.join(self.output_directory, "cantera", file_name) + out_name = os.path.join(self.output_directory, "cantera_from_ck", file_name) if "surface_file" in kwargs: out_name = out_name.replace("-gas.", ".") cantera_dir = os.path.dirname(out_name) diff --git a/rmgpy/species.py b/rmgpy/species.py index e0ae40a2f0..32d8d8d425 100644 --- a/rmgpy/species.py +++ b/rmgpy/species.py @@ -443,9 +443,17 @@ def to_cantera(self, use_chemkin_identifier=False): else: element_dict[symbol] += 1 if use_chemkin_identifier: - ct_species = ct.Species(self.to_chemkin(), element_dict) + label = self.to_chemkin() else: - ct_species = ct.Species(self.label, element_dict) + label = self.label + + if self.contains_surface_site() and element_dict["X"] > 1: + # for multidentate adsorbates, 'size' is the same as 'sites'? for some reason, cantera won't take the input 'size,' so will need to use 'sites' + ct_species = ct.Species(label, element_dict, size=element_dict["X"]) + # hopefully this will be fixed soon, so that ct.Species can take a 'sites' parameter or that cantera can read input files with 'size' specified + else: + ct_species = ct.Species(label, element_dict) + if self.thermo: try: ct_species.thermo = self.thermo.to_cantera() diff --git a/rmgpy/yaml_cantera.py b/rmgpy/yaml_cantera.py new file mode 100644 index 0000000000..c4684bf320 --- /dev/null +++ b/rmgpy/yaml_cantera.py @@ -0,0 +1,378 @@ +#!/usr/bin/env python3 + +############################################################################### +# # +# RMG - Reaction Mechanism Generator # +# # +# Copyright (c) 2002-2024 Prof. William H. Green (whgreen@mit.edu), # +# Prof. Richard H. West (r.west@neu.edu) and the RMG Team (rmg_dev@mit.edu) # +# # +# Permission is hereby granted, free of charge, to any person obtaining a # +# copy of this software and associated documentation files (the 'Software'), # +# to deal in the Software without restriction, including without limitation # +# the rights to use, copy, modify, merge, publish, distribute, sublicense, # +# and/or sell copies of the Software, and to permit persons to whom the # +# Software is furnished to do so, subject to the following conditions: # +# # +# The above copyright notice and this permission notice shall be included in # +# all copies or substantial portions of the Software. # +# # +# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING # +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER # +# DEALINGS IN THE SOFTWARE. # +# # +############################################################################### + +""" +This file defines functions for outputting the RMG generated mechanism to +a yaml file that can be read by Cantera +""" + + +import os +import yaml + +from rmgpy.species import Species +from rmgpy.kinetics.arrhenius import ( + MultiArrhenius, + MultiPDepArrhenius, +) +from rmgpy.util import make_output_subdirectory +from datetime import datetime +from rmgpy.chemkin import get_species_identifier + + +def write_cantera( + spcs, + rxns, + surface_site_density=None, + solvent=None, + solvent_data=None, + path="chem.yml", +): + """ + Writes yaml file depending on the type of system (gas-phase, catalysis). + Writes beginning lines of yaml file, then uses yaml.dump(result_dict) to write species/reactions info. + """ + + # intro to file will change depending on the presence of surface species + is_surface = False + for spc in spcs: + if spc.contains_surface_site(): + is_surface = True + if is_surface: + result_dict = get_mech_dict_surface( + spcs, rxns, solvent=solvent, solvent_data=solvent_data + ) + phases_block, elements_block = get_phases_elements_with_surface( + spcs, surface_site_density + ) + else: + result_dict = get_mech_dict_nonsurface( + spcs, rxns, solvent=solvent, solvent_data=solvent_data + ) + phases_block, elements_block = get_phases_elements_gas_only(spcs) + + with open(path, "w") as f: + # generator line + f.write("generator: RMG\n") + + # datetime object containing current date and time + now = datetime.now() + dt_string = now.strftime("%a, %d %b %Y %H:%M:%S") + f.write(f"date: {dt_string}\n") + + # units line + f.write( + "\nunits: {length: cm, time: s, quantity: mol, activation-energy: kcal/mol}\n\n" + ) + + f.write(phases_block) + f.write(elements_block) + + yaml.dump(result_dict, stream=f, sort_keys=False) + + +def get_phases_elements_gas_only(spcs): + """ + Returns 'phases' and 'elements' sections for a file + with only gas-phase species/reactions. + """ + sorted_species = sorted(spcs, key=lambda spcs: spcs.index) + species_to_write = [get_species_identifier(spec) for spec in sorted_species] + # make sure species with "[" or "]" is in quotes + species_to_write = [ + f"'{s}'" if "[" in s or "{" in s or "]" in s or "}" in s else s + for s in species_to_write + ] + phases_block = f""" +phases: +- name: gas + thermo: ideal-gas + elements: [H, D, T, C, Ci, O, Oi, N, Ne, Ar, He, Si, S, F, Cl, Br, I] + species: [{', '.join(species_to_write)}] + kinetics: gas + transport: mixture-averaged + state: {{T: 300.0, P: 1 atm}} +""" + + elements_block = """ +elements: +- symbol: Ci + atomic-weight: 13.003 +- symbol: D + atomic-weight: 2.014 +- symbol: Oi + atomic-weight: 17.999 +- symbol: T + atomic-weight: 3.016 + +""" + return phases_block, elements_block + + +def get_phases_elements_with_surface(spcs, surface_site_density): + """ + Yaml files with surface species begin with the following blocks of text, + which includes TWO phases instead of just one. + Returns 'phases' and 'elements' sections. + """ + surface_species = [] + gas_species = [] + for spc in spcs: + + if spc.contains_surface_site(): + surface_species.append(spc) + else: + gas_species.append(spc) + + sorted_surface_species = sorted( + surface_species, key=lambda surface_species: surface_species.index + ) + + surface_species_to_write = [ + get_species_identifier(s) for s in sorted_surface_species + ] + + # make sure species with "[" or "]" is in quotes + surface_species_to_write = [ + f"'{s}'" if "[" in s or "{" in s or "]" in s or "}" in s else s + for s in surface_species_to_write + ] + + sorted_gas_species = sorted(gas_species, key=lambda gas_species: gas_species.index) + gas_species_to_write = [get_species_identifier(s) for s in sorted_gas_species] + + # make sure species with "[" or "]" is in quotes + gas_species_to_write = [ + f"'{s}'" if "[" in s or "{" in s or "]" in s or "}" in s else s + for s in gas_species_to_write + ] + + phases_block = f""" +phases: +- name: gas + thermo: ideal-gas + elements: [H, D, T, C, Ci, O, Oi, N, Ne, Ar, He, Si, S, F, Cl, Br, I] + species: [{', '.join(gas_species_to_write)}] + kinetics: gas + reactions: [gas_reactions] + transport: mixture-averaged + state: {{T: 300.0, P: 1 atm}} + +- name: {surface_species[0].smiles.replace("[","").replace("]","")}_surface + thermo: ideal-surface + adjacent-phases: [gas] + elements: [H, D, T, C, Ci, O, Oi, N, Ne, Ar, He, Si, S, F, Cl, Br, I, X] + species: [{', '.join(surface_species_to_write)}] + kinetics: surface + reactions: [surface_reactions] + site-density: {surface_site_density * 1e-4 } +""" + # surface_site_density * 1e-4 #in units of mol/cm^2 + + elements_block = """ +elements: +- symbol: Ci + atomic-weight: 13.003 +- symbol: D + atomic-weight: 2.014 +- symbol: Oi + atomic-weight: 17.999 +- symbol: T + atomic-weight: 3.016 +- symbol: X + atomic-weight: 195.083 + +""" + return phases_block, elements_block + + +def get_mech_dict_surface(spcs, rxns, solvent="solvent", solvent_data=None): + """ + For systems with surface species/reactions. + Adds 'species', 'gas-reactions', and 'surface-reactions' to result_dict. + """ + gas_rxns = [] + surface_rxns = [] + for rxn in rxns: + if rxn.is_surface_reaction(): + surface_rxns.append(rxn) + else: + gas_rxns.append(rxn) + + names = [x.label for x in spcs] + for i, name in enumerate(names): # fix duplicate names + if names.count(name) > 1: + names[i] += "-" + str(names.count(name)) + + result_dict = dict() + result_dict["species"] = [species_to_dict(x, spcs, names=names) for x in spcs] + + # separate gas and surface reactions + + gas_reactions = [] + for rmg_rxn in gas_rxns: + gas_reactions.extend(reaction_to_dicts(rmg_rxn, spcs)) + result_dict["gas_reactions"] = gas_reactions + + surface_reactions = [] + for rmg_rxn in surface_rxns: + surface_reactions.extend(reaction_to_dicts(rmg_rxn, spcs)) + result_dict["surface_reactions"] = surface_reactions + + return result_dict + + +def get_mech_dict_nonsurface(spcs, rxns, solvent="solvent", solvent_data=None): + """ + For gas-phase systems. + Adds 'species' and 'reactions' to result_dict. + """ + names = [x.label for x in spcs] + for i, name in enumerate(names): # fix duplicate names + if names.count(name) > 1: + names[i] += "-" + str(names.count(name)) + + result_dict = dict() + result_dict["species"] = [species_to_dict(x, spcs, names=names) for x in spcs] + + reactions = [] + for rmg_rxn in rxns: + reactions.extend(reaction_to_dicts(rmg_rxn, spcs)) + result_dict["reactions"] = reactions + + return result_dict + + +def reaction_to_dicts(obj, spcs): + """ + Takes an RMG reaction object (obj), returns a list of dictionaries + for YAML properties. For most reaction objects the list will be of + length 1, but a MultiArrhenius or MultiPDepArrhenius will be longer. + """ + + reaction_list = [] + if isinstance(obj.kinetics, MultiArrhenius) or isinstance( + obj.kinetics, MultiPDepArrhenius + ): + list_of_cantera_reactions = obj.to_cantera(use_chemkin_identifier=True) + else: + list_of_cantera_reactions = [obj.to_cantera(use_chemkin_identifier=True)] + + for reaction in list_of_cantera_reactions: + reaction_data = reaction.input_data + efficiencies = getattr(obj.kinetics, "efficiencies", {}) + if efficiencies: + reaction_data["efficiencies"] = { + spcs[i].to_chemkin(): float(val) + for i, val in enumerate( + obj.kinetics.get_effective_collider_efficiencies(spcs) + ) + if val != 1 + } + reaction_list.append(reaction_data) + + return reaction_list + + +def species_to_dict(obj, spc, names=None, label="solvent"): + """ + Takes an RMG species object (obj), returns a list of dictionaries + for YAML properties. Also adds in the number of surface sites + ('sites') to dictionary. + """ + + result_dict = dict() + + if isinstance(obj, Species): + s = obj.to_cantera(use_chemkin_identifier=True) + species_data = s.input_data + try: + result_dict["note"] = obj.transport_data.comment + except: + pass + if "size" in species_data: + sites = species_data["size"] + species_data.pop("size", None) + species_data["sites"] = sites + species_data.update(result_dict) + return ( + species_data # returns composition, name, thermo, and transport, and note + ) + else: + raise Exception("Species object must be an RMG Species object") + + +class CanteraWriter(object): + """ + This class listens to a RMG subject + and writes an YAML file with the current state of the RMG model, + to a yaml subfolder. + + + A new instance of the class can be appended to a subject as follows: + + rmg = ... + listener = CanteraWriter(outputDirectory) + rmg.attach(listener) + + Whenever the subject calls the .notify() method, the + .update() method of the listener will be called. + + To stop listening to the subject, the class can be detached + from its subject: + + rmg.detach(listener) + + """ + + def __init__(self, output_directory=""): + super(CanteraWriter, self).__init__() + self.output_directory = output_directory + make_output_subdirectory(output_directory, "cantera") + + def update(self, rmg): + + solvent_data = None + if rmg.solvent: + solvent_data = rmg.database.solvation.get_solvent_data(rmg.solvent) + + surface_site_density = None + if rmg.reaction_model.surface_site_density: + surface_site_density = rmg.reaction_model.surface_site_density.value_si + + write_cantera( + rmg.reaction_model.core.species, + rmg.reaction_model.core.reactions, + surface_site_density=surface_site_density, + solvent=rmg.solvent, + solvent_data=solvent_data, + path=os.path.join(self.output_directory, "cantera", "chem{}.yaml").format( + len(rmg.reaction_model.core.species) + ), + ) diff --git a/rmgpy/yml.py b/rmgpy/yaml_rms.py similarity index 98% rename from rmgpy/yml.py rename to rmgpy/yaml_rms.py index cd8b9de5b5..f5b1ecc7c1 100644 --- a/rmgpy/yml.py +++ b/rmgpy/yaml_rms.py @@ -48,7 +48,7 @@ from rmgpy.util import make_output_subdirectory -def convert_chemkin_to_yml(chemkin_path, dictionary_path=None, output="chem.rms"): +def convert_chemkin_to_rms(chemkin_path, dictionary_path=None, output="chem.rms"): if dictionary_path: spcs, rxns = load_chemkin_file(chemkin_path, dictionary_path=dictionary_path) else: @@ -56,7 +56,7 @@ def convert_chemkin_to_yml(chemkin_path, dictionary_path=None, output="chem.rms" write_yml(spcs, rxns, path=output) -def write_yml(spcs, rxns, solvent=None, solvent_data=None, path="chem.yml"): +def write_rms(spcs, rxns, solvent=None, solvent_data=None, path="chem.rms"): result_dict = get_mech_dict(spcs, rxns, solvent=solvent, solvent_data=solvent_data) with open(path, 'w') as f: yaml.dump(result_dict, stream=f) @@ -253,5 +253,5 @@ def update(self, rmg): solvent_data = None if rmg.solvent: solvent_data = rmg.database.solvation.get_solvent_data(rmg.solvent) - write_yml(rmg.reaction_model.core.species, rmg.reaction_model.core.reactions, solvent=rmg.solvent, solvent_data=solvent_data, + write_rms(rmg.reaction_model.core.species, rmg.reaction_model.core.reactions, solvent=rmg.solvent, solvent_data=solvent_data, path=os.path.join(self.output_directory, 'rms', 'chem{}.rms').format(len(rmg.reaction_model.core.species))) diff --git a/test/rmgpy/chemkinTest.py b/test/rmgpy/chemkinTest.py index 09099b865c..69692ea4bb 100644 --- a/test/rmgpy/chemkinTest.py +++ b/test/rmgpy/chemkinTest.py @@ -805,9 +805,9 @@ def test_write_bidentate_species(self): chemkin_path = os.path.join(folder, "surface", "chem-surface.inp") dictionary_path = os.path.join(folder, "surface", "species_dictionary.txt") chemkin_save_path = os.path.join(folder, "surface", "chem-surface-test.inp") - species, reactions = load_chemkin_file(chemkin_path, dictionary_path) - - surface_atom_count = species[3].molecule[0].get_num_atoms("X") + species, reactions = load_chemkin_file(chemkin_path, dictionary_path, use_chemkin_names=True) + chox3 = next(iter(s for s in species if s.label=="CHOX3")) + surface_atom_count = chox3.molecule[0].get_num_atoms("X") assert surface_atom_count == 3 save_chemkin_surface_file( chemkin_save_path, @@ -817,17 +817,12 @@ def test_write_bidentate_species(self): check_for_duplicates=False, ) - bidentate_test = " CH2OX2(52)/2/ \n" - tridentate_test = " CHOX3(61)/3/ \n" + bidentate_test = "CH2OX2/2/" + tridentate_test = "CHOX3/3/" with open(chemkin_save_path, "r") as f: - for i, line in enumerate(f): - if i == 3: - bidentate_read = line - if i == 4: - tridentate_read = line - - assert bidentate_test.strip() == bidentate_read.strip() - assert tridentate_test.strip() == tridentate_read.strip() + lines = [line.strip() for line in f] + assert bidentate_test in lines + assert tridentate_test in lines os.remove(chemkin_save_path) diff --git a/test/rmgpy/test_data/chemkin/chemkin_py/surface/chem-gas.inp b/test/rmgpy/test_data/chemkin/chemkin_py/surface/chem-gas.inp new file mode 100644 index 0000000000..0815daee86 --- /dev/null +++ b/test/rmgpy/test_data/chemkin/chemkin_py/surface/chem-gas.inp @@ -0,0 +1,54 @@ +ELEMENTS H C O N Ar END + +SPECIES + Ar + N2 + ethane + CH3 + CH4 + O2 +END + + +THERM ALL + 300.000 1000.000 5000.000 + +Ar Ar1 G200.000 6000.000 1000.00 1 + 2.50000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 2 +-7.45375000E+02 4.37967000E+00 2.50000000E+00 0.00000000E+00 0.00000000E+00 3 + 0.00000000E+00 0.00000000E+00-7.45375000E+02 4.37967000E+00 4 + +N2 N 2 G200.000 6000.000 1000.00 1 + 2.95258000E+00 1.39690000E-03-4.92632000E-07 7.86010000E-11-4.60755000E-15 2 +-9.23949000E+02 5.87189000E+00 3.53101000E+00-1.23661000E-04-5.02999000E-07 3 + 2.43531000E-09-1.40881000E-12-1.04698000E+03 2.96747000E+00 4 + +ethane C 2 H 6 G100.000 5000.000 954.51 1 + 4.58979542E+00 1.41508364E-02-4.75965787E-06 8.60302924E-10-6.21723861E-14 2 +-1.27217507E+04-3.61718979E+00 3.78034578E+00-3.24276131E-03 5.52385395E-05 3 +-6.38587729E-08 2.28639990E-11-1.16203414E+04 5.21029728E+00 4 + +CH3 C 1 H 3 G100.000 5000.000 1337.62 1 + 3.54143640E+00 4.76790043E-03-1.82150130E-06 3.28880372E-10-2.22548587E-14 2 + 1.62239681E+04 1.66047100E+00 3.91546905E+00 1.84152818E-03 3.48746163E-06 3 +-3.32752224E-09 8.49972445E-13 1.62856393E+04 3.51736167E-01 4 + +CH4 C 1H 4 G 100.000 5000.000 1084.12 1 + 9.08256809E-01 1.14541005E-02-4.57174656E-06 8.29193594E-10-5.66316470E-14 2 +-9.71997053E+03 1.39931449E+01 4.20541679E+00-5.35559144E-03 2.51123865E-05 3 +-2.13763581E-08 5.97526898E-12-1.01619434E+04-9.21284857E-01 4 + +O2 O 2 G 100.000 5000.000 1074.56 1 + 3.15382774E+00 1.67803224E-03-7.69967755E-07 1.51273954E-10-1.08781177E-14 2 +-1.04082031E+03 6.16751905E+00 3.53732118E+00-1.21570202E-03 5.31615358E-06 3 +-4.89440364E-09 1.45843807E-12-1.03858843E+03 4.68368633E+00 4 + +END + + +REACTIONS KCAL/MOLE MOLES + +CH3 +CH3 = ethane 8.260e+17 -1.400 1.000 + +END + diff --git a/test/rmgpy/test_data/chemkin/chemkin_py/surface/chem-surface.inp b/test/rmgpy/test_data/chemkin/chemkin_py/surface/chem-surface.inp index 29dfedd651..992e28bb70 100644 --- a/test/rmgpy/test_data/chemkin/chemkin_py/surface/chem-surface.inp +++ b/test/rmgpy/test_data/chemkin/chemkin_py/surface/chem-surface.inp @@ -1,61 +1,63 @@ ELEMENTS H - D /2.014/ - T /3.016/ C - CI /13.003/ O - OI /18.000/ - N - Ne Ar - He - Si - S - F - Cl - Br - I X /195.083/ END -SPECIES - X(1) - H*(10) - CH2OX2(52) - CHOX3(61) +SITE SDEN/2.4830E-09/ ! mol/cm^2 + X + HX + OX + CH4X + CH2OX2/2/ + CHOX3/3/ END THERM ALL 300.000 1000.000 5000.000 -X(1) X 1 G 100.000 5000.000 1554.80 1 +X X 1 G 100.000 5000.000 1554.80 1 1.60299900E-01-2.52235409E-04 1.14181275E-07-1.21471653E-11 3.85790025E-16 2 -7.08100885E+01-9.09527530E-01 7.10139498E-03-4.25619522E-05 8.98533016E-08 3 -7.80193649E-11 2.32465471E-14-8.76101712E-01-3.11211229E-02 4 -H*(10) H 1X 1 G 100.000 5000.000 952.91 1 +HX H 1X 1 G 100.000 5000.000 952.91 1 2.80339655E+00-5.41047017E-04 4.99507978E-07-7.54963647E-11 3.06772366E-15 2 -2.34636021E+03-1.59436787E+01-3.80965452E-01 5.47228709E-03 2.60912778E-06 3 -9.64961980E-09 4.63946753E-12-1.40561079E+03 1.01725550E+00 4 -CH2OX2(52) C 1H 2O 1X 2G 100.000 5000.000 777.56 1 +CH2OX2 C 1H 2O 1X 2G 100.000 5000.000 777.56 1 3.18259452E+00 1.04657053E-02-5.00464545E-06 9.79545329E-10-6.90993489E-14 2 -1.82783093E+04-1.50760119E+01-1.44645549E+00 3.42803139E-02-5.09483795E-05 3 4.03732306E-08-1.27356452E-11-1.75584787E+04 6.09156343E+00 4 -CHOX3(61) C 1H 1O 1X 3G 100.000 5000.000 689.33 1 +CHOX3 C 1H 1O 1X 3G 100.000 5000.000 689.33 1 3.47343833E+00 7.03348332E-03-3.51779073E-06 6.98047945E-10-4.93382204E-14 2 -1.53390119E+04-1.79870016E+01-1.32716780E+00 3.17587333E-02-5.05065871E-05 3 3.95520303E-08-1.17505741E-11-1.46027731E+04 3.92679635E+00 4 -END +OX O 1X 1 G 298.000 2000.000 1000.00 1 + 2.90244691E+00-3.38584457E-04 6.43372619E-07-3.66326660E-10 6.90093884E-14 2 +-1.70497471E+04-1.52559728E+01-2.94475701E-01 1.44162624E-02-2.61322704E-05 3 + 2.19005957E-08-6.98019420E-12-1.64619234E+04-1.99445623E-01 4 +CH4X C 1H 4X 1 G 298.000 2000.000 1000.00 1 + 9.54139378E+00-1.04025134E-02 1.83777400E-05-9.66765130E-09 1.71211379E-12 2 +-1.34475614E+04-3.55638434E+01 4.85496247E+00-5.54134984E-03 3.01198105E-05 3 +-2.99225917E-08 1.00502514E-11-1.17096278E+04-9.25620913E+00 4 +END REACTIONS KCAL/MOLE MOLES -X(1)+X(1)+CH2OX2(52)=H*(10)+CHOX3(61) 7.420000e+21 0.000 3.058 +X + X + CH2OX2 = HX + CHOX3 7.420000e+21 0.000 3.058 + +X + CH4 <=> CH4X 8.000e-03 0.000 0.000 + STICK + +OX + OX <=> X + X + O2 3.700000e+21 0.000 66.611 END diff --git a/test/rmgpy/test_data/chemkin/chemkin_py/surface/species_dictionary.txt b/test/rmgpy/test_data/chemkin/chemkin_py/surface/species_dictionary.txt index f6acc21000..c03ebd7541 100644 --- a/test/rmgpy/test_data/chemkin/chemkin_py/surface/species_dictionary.txt +++ b/test/rmgpy/test_data/chemkin/chemkin_py/surface/species_dictionary.txt @@ -1,11 +1,11 @@ -X(1) +X 1 X u0 p0 c0 -H*(10) +HX 1 H u0 p0 c0 {2,S} 2 X u0 p0 c0 {1,S} -CHOX3(61) +CHOX3 1 O u0 p2 c0 {2,S} {5,S} 2 C u0 p0 c0 {1,S} {3,S} {4,S} {6,S} 3 H u0 p0 c0 {2,S} @@ -13,10 +13,58 @@ CHOX3(61) 5 X u0 p0 c0 {1,S} 6 X u0 p0 c0 {2,S} -CH2OX2(52) +CH2OX2 1 O u0 p2 c0 {2,S} {6,S} 2 C u0 p0 c0 {1,S} {3,S} {4,S} {5,S} 3 H u0 p0 c0 {2,S} 4 H u0 p0 c0 {2,S} 5 X u0 p0 c0 {2,S} -6 X u0 p0 c0 {1,S} \ No newline at end of file +6 X u0 p0 c0 {1,S} + +Ar +1 Ar u0 p4 c0 + +N2 +1 N u0 p1 c0 {2,T} +2 N u0 p1 c0 {1,T} + +ethane +1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S} +2 C u0 p0 c0 {1,S} {6,S} {7,S} {8,S} +3 H u0 p0 c0 {1,S} +4 H u0 p0 c0 {1,S} +5 H u0 p0 c0 {1,S} +6 H u0 p0 c0 {2,S} +7 H u0 p0 c0 {2,S} +8 H u0 p0 c0 {2,S} + +CH3 +multiplicity 2 +1 C u1 p0 c0 {2,S} {3,S} {4,S} +2 H u0 p0 c0 {1,S} +3 H u0 p0 c0 {1,S} +4 H u0 p0 c0 {1,S} + +O2 +multiplicity 3 +1 O u1 p2 c0 {2,S} +2 O u1 p2 c0 {1,S} + +CH4 +1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S} +2 H u0 p0 c0 {1,S} +3 H u0 p0 c0 {1,S} +4 H u0 p0 c0 {1,S} +5 H u0 p0 c0 {1,S} + +CH4X +1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S} +2 H u0 p0 c0 {1,S} +3 H u0 p0 c0 {1,S} +4 H u0 p0 c0 {1,S} +5 H u0 p0 c0 {1,S} +6 X u0 p0 c0 + +OX +1 O u0 p2 c0 {2,D} +2 X u0 p0 c0 {1,D} \ No newline at end of file diff --git a/test/rmgpy/tools/canteramodelTest.py b/test/rmgpy/tools/canteramodelTest.py index c2e1b1d57c..81b7d6844e 100644 --- a/test/rmgpy/tools/canteramodelTest.py +++ b/test/rmgpy/tools/canteramodelTest.py @@ -104,6 +104,27 @@ def setup_class(self): self.ctSpecies = job.model.species() self.ctReactions = job.model.reactions() + # Now load surface species and kinetics + folder = os.path.join(os.path.dirname(os.path.dirname(__file__)), "test_data", "chemkin", "chemkin_py") + chemkin_path = os.path.join(folder, "surface", "chem-gas.inp") + chemkin_surface_path = os.path.join(folder, "surface", "chem-surface.inp") + dictionary_path = os.path.join(folder, "surface", "species_dictionary.txt") + species, reactions = load_chemkin_file(chemkin_surface_path, dictionary_path) + self.rmg_surface_ct_species = [spec.to_cantera(use_chemkin_identifier=True) for spec in species] + self.rmg_surface_ct_reactions = [] + for rxn in reactions: + converted_reactions = rxn.to_cantera(species, use_chemkin_identifier=True) + if isinstance(converted_reactions, list): + self.rmg_surface_ct_reactions.extend(converted_reactions) + else: + self.rmg_surface_ct_reactions.append(converted_reactions) + job = Cantera() + job.surface = True + job.load_chemkin_model(chemkin_path, surface_file=chemkin_surface_path, quiet=True) + self.ct_surface_species = job.surface.species() + self.ct_surface_reactions = job.surface.reactions() + + def test_species_conversion(self): """ Test that species objects convert properly @@ -115,9 +136,31 @@ def test_species_conversion(self): def test_reaction_conversion(self): """ - Test that species objects convert properly + Test that reaction objects convert properly """ from rmgpy.tools.canteramodel import check_equivalent_cantera_reaction for i in range(len(self.ctReactions)): assert check_equivalent_cantera_reaction(self.ctReactions[i], self.rmg_ctReactions[i]) + + def test_surface_species_conversion(self): + """ + Test that surface species objects convert properly + """ + from rmgpy.tools.canteramodel import check_equivalent_cantera_species + + for i in range(len(self.ct_surface_species)): + #print("Chemkin-to-Cantera:", self.ct_surfaceSpecies[i].input_data) + #print("Chemkin-to-RMG-to-Cantera:", self.rmg_surface_ctSpecies[i].input_data) + assert check_equivalent_cantera_species(self.ct_surface_species[i], self.rmg_surface_ct_species[i]) + + def test_surface_reaction_conversion(self): + """ + Test that surface reaction objects convert properly + """ + from rmgpy.tools.canteramodel import check_equivalent_cantera_reaction + + for i in range(len(self.ct_surface_reactions)): + #print("Chemkin-to-Cantera:", self.ct_surfaceReactions[i].input_data) + #print("Chemkin-to-RMG-to-Cantera:", self.rmg_surface_ctReactions[i].input_data) + assert check_equivalent_cantera_reaction(self.ct_surface_reactions[i], self.rmg_surface_ct_reactions[i]) \ No newline at end of file