Create YAML file outputs in RMG for use in Cantera

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)
 
 

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()
 
 ################################################################################
 

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):

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 = {}

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)

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()