Add and improve to_cantera methods for conversion to in-memory Cantera objects

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

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

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)
 

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
+