Yaml writer addition

rmgpy/kinetics/arrhenius.pyx

@@ -946,17 +946,24 @@ cdef class PDepArrhenius(PDepKineticsModel):
     def set_cantera_kinetics(self, ct_reaction, species_list):
         """
         Sets a Cantera PlogReaction()'s `rates` attribute with
-        A list of tuples containing [(pressure in Pa, cantera arrhenius object), (..)]
+        a list of tuples containing [(pressure in Pa, cantera arrhenius object), ...].
+
+        A ``MultiArrhenius`` entry (from chemkin PLOG blocks with duplicate
+        pressures) is expanded into one tuple per inner Arrhenius; Cantera's
+        PlogRate sums duplicate-pressure entries at evaluation.
         """
         import cantera as ct
-        import copy
         assert isinstance(ct_reaction.rate, ct.PlogRate), "Must have a Cantera PlogRate attribute"
 
-        pressures = copy.deepcopy(self._pressures.value_si)
-        ctArrhenius = [arr.to_cantera_kinetics(arrhenius_class=True) for arr in self.arrhenius]
+        rate_pairs = []
+        for P, arr in zip(self._pressures.value_si, self.arrhenius):
+            if isinstance(arr, MultiArrhenius):
+                for sub in arr.arrhenius:
+                    rate_pairs.append((P, sub.to_cantera_kinetics(arrhenius_class=True)))
+            else:
+                rate_pairs.append((P, arr.to_cantera_kinetics(arrhenius_class=True)))
 
-        new_rates = ct.PlogRate(list(zip(pressures, ctArrhenius)))
-        ct_reaction.rate = new_rates
+        ct_reaction.rate = ct.PlogRate(rate_pairs)
 
 ################################################################################
 

rmgpy/reaction.py

@@ -425,6 +425,19 @@ def to_cantera(self, species_list=None, use_chemkin_identifier=False):
         # Now we set the kinetics.
         self.kinetics.set_cantera_kinetics(ct_reaction, species_list)
 
+        # Coverage dependencies are not handled by set_cantera_kinetics; set them here.
+        # Cantera's coverage_dependencies E is in J/kmol; RMG's value_si is J/mol.
+        if hasattr(self.kinetics, 'coverage_dependence') and self.kinetics.coverage_dependence:
+            cov_deps = {}
+            for sp, params in self.kinetics.coverage_dependence.items():
+                sp_label = sp.to_chemkin() if use_chemkin_identifier else sp.label
+                cov_deps[sp_label] = {
+                    'a': params['a'].value_si,
+                    'm': params['m'].value_si,
+                    'E': params['E'].value_si * 1000,  # J/mol → J/kmol
+                }
+            ct_reaction.rate.coverage_dependencies = cov_deps
+
         return ct_reaction
 
     def get_url(self):

rmgpy/rmg/main.py

@@ -81,8 +81,11 @@
 from rmgpy.stats import ExecutionStatsWriter
 from rmgpy.thermo.thermoengine import submit
 from rmgpy.tools.plot import plot_sensitivity
+from rmgpy.tools.compare_cantera_yaml import compare_yaml_files, compare_yaml_files_and_report
 from rmgpy.tools.uncertainty import Uncertainty, process_local_results
-from rmgpy.yml import RMSWriter
+from rmgpy.yaml_rms import RMSWriter
+from rmgpy.yaml_cantera1 import CanteraWriter1
+from rmgpy.yaml_cantera2 import CanteraWriter2
 
 ################################################################################
 
@@ -784,7 +787,8 @@ def register_listeners(self, requires_rms=False):
         self.attach(ChemkinWriter(self.output_directory))
 
         self.attach(RMSWriter(self.output_directory))
-
+        self.attach(CanteraWriter1(self.output_directory))
+        self.attach(CanteraWriter2(self.output_directory))
         if self.generate_output_html:
             self.attach(OutputHTMLWriter(self.output_directory))
 
@@ -1232,15 +1236,17 @@ def execute(self, initialize=True, **kwargs):
 
         self.run_model_analysis()
 
-        # generate Cantera files chem.yaml & chem_annotated.yaml in a designated `cantera` output folder
+        # generate Cantera files chem.yaml & chem_annotated.yaml in designated Cantera output folders
         try:
+            logging.info("Translating final chemkin file into Cantera yaml.")
+            translated_cantera_file = None
             if any([s.contains_surface_site() for s in self.reaction_model.core.species]):
                 # Surface (catalytic) chemistry
-                self.generate_cantera_files(
+                translated_cantera_file = self.generate_cantera_files_from_chemkin(
                     os.path.join(self.output_directory, "chemkin", "chem-gas.inp"),
                     surface_file=(os.path.join(self.output_directory, "chemkin", "chem-surface.inp")),
                 )
-                self.generate_cantera_files(
+                self.generate_cantera_files_from_chemkin(
                     os.path.join(self.output_directory, "chemkin", "chem_annotated-gas.inp"),
                     surface_file=(os.path.join(self.output_directory, "chemkin", "chem_annotated-surface.inp")),
                 )
@@ -1273,8 +1279,22 @@ def execute(self, initialize=True, **kwargs):
                         _add_coverage_dependence_to_cantera_yaml(yaml_path, coverage_deps)
 
             else:  # gas phase only
-                self.generate_cantera_files(os.path.join(self.output_directory, "chemkin", "chem.inp"))
-                self.generate_cantera_files(os.path.join(self.output_directory, "chemkin", "chem_annotated.inp"))
+                translated_cantera_file = self.generate_cantera_files_from_chemkin(
+                    os.path.join(self.output_directory, "chemkin", "chem.inp")
+                )
+                self.generate_cantera_files_from_chemkin(
+                    os.path.join(self.output_directory, "chemkin", "chem_annotated.inp")
+                )
+
+            # Compare translated Cantera files and directly generated Cantera files
+
+            compare_yaml_files_and_report(translated_cantera_file,
+                                          os.path.join(self.output_directory, "cantera1", "chem.yaml"),
+                                          output=os.path.join(self.output_directory, "cantera1", "comparison_report.txt"))
+            compare_yaml_files_and_report(translated_cantera_file,
+                                          os.path.join(self.output_directory, "cantera2", "chem.yaml"),
+                                          output=os.path.join(self.output_directory, "cantera2", "comparison_report.txt"))
+
         except EnvironmentError:
             logging.exception("Could not generate Cantera files due to EnvironmentError. Check read\\write privileges in output directory.")
         except Exception:
@@ -1845,14 +1865,15 @@ def process_reactions_to_species(self, obj):
             raise TypeError("improper call, obj input was incorrect")
         return potential_spcs
 
-    def generate_cantera_files(self, chemkin_file, **kwargs):
+    def generate_cantera_files_from_chemkin(self, chemkin_file, **kwargs):
         """
         Convert a chemkin mechanism chem.inp file to a cantera mechanism file chem.yaml
-        and save it in the cantera directory
+        and save it in the cantera directory. 
+        Returns the path to the generated cantera file.
         """
         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)
@@ -1870,6 +1891,7 @@ def generate_cantera_files(self, chemkin_file, **kwargs):
             logging.exception("Error converting to Cantera format.")
             logging.info("Trying again without transport data file.")
             parser.convert_mech(chemkin_file, out_name=out_name, quiet=True, permissive=True, **kwargs)
+        return out_name
 
     def initialize_reaction_threshold_and_react_flags(self):
         num_core_species = len(self.reaction_model.core.species)
@@ -2157,8 +2179,9 @@ def get_git_commit(self, module_path):
 
         if os.path.exists(os.path.join(module_path, "..", ".git")):
             try:
-                return subprocess.check_output(["git", "log", "--format=%H%n%cd", "-1"], cwd=module_path).splitlines()
-            except:
+                head, date = subprocess.check_output(["git", "log", "--format=%H%n%cd", "-1"], cwd=module_path).splitlines()
+                return head.decode(), date.decode()
+            except (subprocess.CalledProcessError, OSError):
                 return "", ""
         else:
             return "", ""

rmgpy/solver/surface.pyx

@@ -36,7 +36,7 @@ import logging
 cimport cython
 import numpy as np
 cimport numpy as np
-from libc.math cimport exp
+from libc.math cimport exp, pow
 
 import rmgpy.constants as constants
 cimport rmgpy.constants as constants
@@ -478,9 +478,11 @@ cdef class SurfaceReactor(ReactionSystem):
             #: surface species are in mol/m2, gas phase are in mol/m3
             core_species_concentrations[j] = C[j]
 
+        if self.thermo_coverage_dependence or self.coverage_dependence:
+            coverages = np.where(species_on_surface, np.maximum(N / total_sites, 0.0), 0.0)
+
         # Thermodynamic coverage dependence
         if self.thermo_coverage_dependence:
-            coverages = np.where(species_on_surface, N / total_sites, 0.0)
             coverages_squared = coverages * coverages
             temperature_scaled_coverages = -self.T.value_si * coverages
             thermo_dep_coverage = np.empty((6, coverages.shape[0]), dtype=np.float64)
@@ -508,12 +510,12 @@ cdef class SurfaceReactor(ReactionSystem):
             """
             for i, list_of_coverage_deps in self.coverage_dependencies.items():
                 # Species i, Reaction j
-                surface_site_fraction = N[i] / total_sites
-                if surface_site_fraction < 1e-15:
+                surface_site_fraction = coverages[i]
+                if surface_site_fraction <= 1e-6:
                     continue
                 for j, a, m, E in list_of_coverage_deps:
                     coverage_corrections[j] *= 10. ** (a * surface_site_fraction) *\
-                                                surface_site_fraction ** m *\
+                                                pow(surface_site_fraction, m) *\
                                                 exp(-1 * E * surface_site_fraction / (constants.R * self.T.value_si))
             kf = kf * coverage_corrections # make a corrected copy kf, but leave the original array at self.kf unchanged
             kr = kr * coverage_corrections

rmgpy/species.py

@@ -428,13 +428,18 @@ def to_chemkin(self):
         from rmgpy.chemkin import get_species_identifier
         return get_species_identifier(self)
 
-    def to_cantera(self, use_chemkin_identifier=False):
+    def to_cantera(self, use_chemkin_identifier=False, all_species=None):
         """
         Converts the RMG Species object to a Cantera Species object
         with the appropriate thermo data.
 
         If use_chemkin_identifier is set to False, the species label is used
         instead. Be sure that species' labels are unique when setting it False.
+
+        If all_species is provided and this is a surface species with
+        coverage-dependent thermo, the coverage-dependencies are attached to
+        the Cantera Species object via update_user_data() so that they appear
+        in input_data and in Solution.write_yaml() output.
         """
         import cantera as ct
 
@@ -475,6 +480,32 @@ def to_cantera(self, use_chemkin_identifier=False):
         if self.transport_data:
             ct_species.transport = self.transport_data.to_cantera()
 
+        # Attach coverage-dependent thermo if present.
+        # thermo_coverage_dependence keys are adjacency-list strings; we resolve
+        # them to species names here using all_species.  The data is stored via
+        # update_user_data() so it appears in ct_species.input_data and in any
+        # YAML serialised by Cantera (e.g. Solution.write_yaml()).
+        if (all_species is not None
+                and self.contains_surface_site()
+                and self.thermo
+                and hasattr(self.thermo, 'thermo_coverage_dependence')
+                and self.thermo.thermo_coverage_dependence):
+            from rmgpy.molecule.molecule import Molecule
+            cov_deps = {}
+            for adj_list, parameters in self.thermo.thermo_coverage_dependence.items():
+                mol = Molecule().from_adjacency_list(adj_list)
+                for sp in all_species:
+                    if sp.is_isomorphic(mol, strict=False):
+                        dep_name = sp.to_chemkin() if use_chemkin_identifier else sp.label
+                        cov_deps[dep_name] = {
+                            'units': {'energy': 'J', 'quantity': 'mol'},
+                            'enthalpy-coefficients': [v.value_si for v in parameters['enthalpy-coefficients']],
+                            'entropy-coefficients': [v.value_si for v in parameters['entropy-coefficients']],
+                        }
+                        break
+            if cov_deps:
+                ct_species.update_user_data({'coverage-dependencies': cov_deps})
+
         return ct_species
 
     def has_statmech(self):

rmgpy/thermo/nasa.pyx

@@ -421,7 +421,6 @@ cdef class NASA(HeatCapacityModel):
             poly.change_base_entropy(deltaS)
         return self
 
-    # need to modify this to include the thermo coverage dependence
     def to_cantera(self):
         """
         Return the cantera equivalent NasaPoly2 object from this NASA object.