Coverage-dependent thermochemistry for catalysis

rmgpy/rmg/input.py

@@ -120,7 +120,8 @@ def database(
 def catalyst_properties(bindingEnergies=None,
                         surfaceSiteDensity=None,
                         metal=None,
-                        coverageDependence=False):
+                        coverageDependence=False,
+                        thermoCoverageDependence=False,):
     """
     Specify the properties of the catalyst.
     Binding energies of C,H,O,N atoms, and the surface site density.
@@ -167,6 +168,7 @@ def catalyst_properties(bindingEnergies=None,
     else:
         logging.info("Coverage dependence is turned OFF")
     rmg.coverage_dependence = coverageDependence
+    rmg.thermo_coverage_dependence = thermoCoverageDependence
 
 def convert_binding_energies(binding_energies):
     """
@@ -1088,7 +1090,8 @@ def surface_reactor(temperature,
                             sensitive_species=sensitive_species,
                             sensitivity_threshold=sensitivityThreshold,
                             sens_conditions=sens_conditions,
-                            coverage_dependence=rmg.coverage_dependence)
+                            coverage_dependence=rmg.coverage_dependence,
+                            thermo_coverage_dependence=rmg.thermo_coverage_dependence)
     rmg.reaction_systems.append(system)
     system.log_initial_conditions(number=len(rmg.reaction_systems))
 

rmgpy/rmg/main.py

@@ -37,6 +37,7 @@
 import logging
 import marshal
 import os
+import re
 import resource
 import shutil
 import sys
@@ -193,6 +194,7 @@ def clear(self):
         self.surface_site_density = None
         self.binding_energies = None
         self.coverage_dependence = False
+        self.thermo_coverage_dependence = False
         self.forbidden_structures = []
 
         self.reaction_model = None
@@ -510,7 +512,7 @@ def initialize(self, **kwargs):
 
         # Read input file
         self.load_input(self.input_file)
-        
+
         # Check if ReactionMechanismSimulator reactors are being used
         # if RMS is not installed but the user attempted to use it, the load_input_file would have failed
         # if RMS is not installed and they did not use it, we avoid calling certain functions that would raise an error
@@ -523,6 +525,7 @@ def initialize(self, **kwargs):
                 self.reaction_model.core.phase_system.phases["Surface"].site_density = self.surface_site_density.value_si
                 self.reaction_model.edge.phase_system.phases["Surface"].site_density = self.surface_site_density.value_si
         self.reaction_model.coverage_dependence = self.coverage_dependence
+        self.reaction_model.thermo_coverage_dependence = self.thermo_coverage_dependence
 
         if kwargs.get("restart", ""):
             import rmgpy.rmg.input
@@ -768,7 +771,7 @@ def register_listeners(self, requires_rms=False):
         """
 
         self.attach(ChemkinWriter(self.output_directory))
-        
+
         self.attach(RMSWriter(self.output_directory))
 
         if self.generate_output_html:
@@ -1221,6 +1224,7 @@ def execute(self, initialize=True, **kwargs):
         # generate Cantera files chem.yaml & chem_annotated.yaml in a designated `cantera` output folder
         try:
             if any([s.contains_surface_site() for s in self.reaction_model.core.species]):
+                # Surface (catalytic) chemistry
                 self.generate_cantera_files(
                     os.path.join(self.output_directory, "chemkin", "chem-gas.inp"),
                     surface_file=(os.path.join(self.output_directory, "chemkin", "chem-surface.inp")),
@@ -1229,6 +1233,34 @@ def execute(self, initialize=True, **kwargs):
                     os.path.join(self.output_directory, "chemkin", "chem_annotated-gas.inp"),
                     surface_file=(os.path.join(self.output_directory, "chemkin", "chem_annotated-surface.inp")),
                 )
+
+                if self.thermo_coverage_dependence:
+                    # Build coverage_deps: {species_name: string_to_add_to_yaml}
+                    coverage_deps = {}
+                    for s in self.reaction_model.core.species:
+                        if s.contains_surface_site() and s.thermo.thermo_coverage_dependence:
+                            s_name = s.to_chemkin()
+                            for dep_sp_adj, parameters in s.thermo.thermo_coverage_dependence.items():
+                                mol = Molecule().from_adjacency_list(dep_sp_adj)
+                                for sp in self.reaction_model.core.species:
+                                    if sp.is_isomorphic(mol, strict=False):
+                                        if s_name not in coverage_deps:
+                                            coverage_deps[s_name] = '  coverage-dependencies:'
+                                        coverage_deps[s_name] += f"""
+    {sp.to_chemkin()}:
+      model: {parameters['model']}
+      enthalpy-coefficients: {[v.value_si for v in parameters['enthalpy-coefficients']]}
+      entropy-coefficients: {[v.value_si for v in parameters['entropy-coefficients']]}
+      units: {{energy: J, quantity: mol}}
+"""
+                                        break
+
+                    for yaml_path in [
+                        os.path.join(self.output_directory, "cantera", "chem.yaml"),
+                        os.path.join(self.output_directory, "cantera", "chem_annotated.yaml"),
+                    ]:
+                        _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"))
@@ -2392,6 +2424,45 @@ def obj(y):
             self.scaled_condition_list.append(scaled_new_cond)
         return
 
+def _add_coverage_dependence_to_cantera_yaml(yaml_path, coverage_deps):
+    """Modify a Cantera YAML file in-place to add coverage-dependent surface thermo.
+
+    Makes targeted text insertions rather than loading and re-dumping the whole
+    file, so original formatting is preserved everywhere except the new lines.
+
+    Args:
+        yaml_path: path to the Cantera YAML file to modify
+        coverage_deps: dict mapping species ChemKin names to their coverage-dependency string.
+    """
+    with open(yaml_path, 'r') as f:
+        content = f.read()
+
+    # --- Modify the surface phase ---
+    # Replace 'ideal-surface' with 'coverage-dependent-surface' and add reference-state-coverage.
+    content = content.replace(
+        '  thermo: ideal-surface\n',
+        '  thermo: coverage-dependent-surface\n  reference-state-coverage: 0.11\n',
+        1,
+    )
+
+    # --- Insert coverage-dependencies block after each relevant species entry ---
+    for species_name, deps in coverage_deps.items():
+        match = re.search(r'^- name: ' + re.escape(species_name) + r'\n', content, re.MULTILINE)
+        if not match:
+            logging.warning(
+                f"Species {species_name} not found in {yaml_path}; skipping coverage-dependency insertion."
+            )
+            continue
+
+        after = match.end()
+        end_match = re.search(r'\n(?=(?:- |\n|\w))', content[after:])
+        if end_match:
+            insert_pos = after + end_match.start() + 1
+            content = content[:insert_pos] + deps + content[insert_pos:]
+        else:
+            content = content.rstrip('\n') + '\n' + deps
+    with open(yaml_path, 'w') as f:
+        f.write(content)
 
 def log_conditions(rmg_memories, index):
     """

rmgpy/solver/surface.pyx

@@ -43,6 +43,8 @@ cimport rmgpy.constants as constants
 from rmgpy.quantity import Quantity
 from rmgpy.quantity cimport ScalarQuantity
 from rmgpy.solver.base cimport ReactionSystem
+import copy
+from rmgpy.molecule import Molecule
 
 cdef class SurfaceReactor(ReactionSystem):
     """
@@ -66,9 +68,13 @@ cdef class SurfaceReactor(ReactionSystem):
     cdef public ScalarQuantity surface_site_density
     cdef public np.ndarray reactions_on_surface  # (catalyst surface, not core/edge surface)
     cdef public np.ndarray species_on_surface  # (catalyst surface, not core/edge surface)
+    cdef public np.ndarray thermo_coeff_matrix
+    cdef public np.ndarray stoi_matrix
 
     cdef public bint coverage_dependence
     cdef public dict coverage_dependencies
+    cdef public bint thermo_coverage_dependence
+
 
 
     def __init__(self,
@@ -84,6 +90,7 @@ cdef class SurfaceReactor(ReactionSystem):
                  sensitivity_threshold=1e-3,
                  sens_conditions=None,
                  coverage_dependence=False,
+                 thermo_coverage_dependence=False,
                  ):
         ReactionSystem.__init__(self,
                                 termination,
@@ -103,6 +110,7 @@ cdef class SurfaceReactor(ReactionSystem):
         self.surface_volume_ratio = Quantity(surface_volume_ratio)
         self.surface_site_density = Quantity(surface_site_density)
         self.coverage_dependence = coverage_dependence
+        self.thermo_coverage_dependence = thermo_coverage_dependence
         self.V = 0  # will be set from ideal gas law in initialize_model
         self.constant_volume = True
         self.sens_conditions = sens_conditions
@@ -166,6 +174,10 @@ cdef class SurfaceReactor(ReactionSystem):
                                        )
         cdef np.ndarray[np.int_t, ndim=1] species_on_surface, reactions_on_surface
         cdef Py_ssize_t index
+        cdef np.ndarray thermo_coeff_matrix = np.zeros((len(self.species_index), len(self.species_index), 6), dtype=np.float64)
+        cdef np.ndarray stoi_matrix = np.zeros((self.reactant_indices.shape[0], len(self.species_index)), dtype=np.float64)
+        if self.thermo_coverage_dependence:
+            self.thermo_coeff_matrix = thermo_coeff_matrix
         #: 1 if it's on a surface, 0 if it's in the gas phase
         reactions_on_surface = np.zeros((self.num_core_reactions + self.num_edge_reactions), int)
         species_on_surface = np.zeros((self.num_core_species), int)
@@ -195,6 +207,49 @@ cdef class SurfaceReactor(ReactionSystem):
                 means that Species with index 2 in the current simulation is used in
                 Reaction 3 with parameters a=0.1, m=-1, E=12 kJ/mol
                 """
+        for sp, sp_index in self.species_index.items():
+            if sp.contains_surface_site():
+                if self.thermo_coverage_dependence and sp.thermo.thermo_coverage_dependence:
+                    for spec, parameters in sp.thermo.thermo_coverage_dependence.items():
+                        molecule = Molecule().from_adjacency_list(spec)
+                        for species in self.species_index.keys():
+                            if species.is_isomorphic(molecule, strict=False):
+                                species_index = self.species_index[species]
+                                enthalpy_coeff = np.array([p.value_si for p in parameters['enthalpy-coefficients']])
+                                entropy_coeff = np.array([p.value_si for p in parameters['entropy-coefficients']])
+                                thermo_polynomials = np.concatenate((enthalpy_coeff, entropy_coeff), axis=0)
+                                self.thermo_coeff_matrix[sp_index, species_index] = [x for x in thermo_polynomials]
+        # create a stoichiometry matrix for reaction enthalpy and entropy correction 
+        # due to thermodynamic coverage dependence
+        if self.thermo_coverage_dependence:
+            ir = self.reactant_indices
+            ip = self.product_indices
+            for rxn_id, rxn_stoi_num in enumerate(stoi_matrix):
+                if ir[rxn_id, 0] >= self.num_core_species or ir[rxn_id, 1] >= self.num_core_species or ir[rxn_id, 2] >= self.num_core_species:
+                    continue
+                elif ip[rxn_id, 0] >= self.num_core_species or ip[rxn_id, 1] >= self.num_core_species or ip[rxn_id, 2] >= self.num_core_species:
+                    continue
+                else:
+                    if ir[rxn_id, 1] == -1:  # only one reactant
+                        rxn_stoi_num[ir[rxn_id, 0]] += -1
+                    elif ir[rxn_id, 2] == -1:  # only two reactants
+                        rxn_stoi_num[ir[rxn_id, 0]] += -1
+                        rxn_stoi_num[ir[rxn_id, 1]] += -1
+                    else:  # three reactants
+                        rxn_stoi_num[ir[rxn_id, 0]] += -1
+                        rxn_stoi_num[ir[rxn_id, 1]] += -1
+                        rxn_stoi_num[ir[rxn_id, 2]] += -1
+                    if ip[rxn_id, 1] == -1:  # only one product
+                        rxn_stoi_num[ip[rxn_id, 0]] += 1
+                    elif ip[rxn_id, 2] == -1:  # only two products
+                        rxn_stoi_num[ip[rxn_id, 0]] += 1
+                        rxn_stoi_num[ip[rxn_id, 1]] += 1
+                    else:  # three products
+                        rxn_stoi_num[ip[rxn_id, 0]] += 1
+                        rxn_stoi_num[ip[rxn_id, 1]] += 1
+                        rxn_stoi_num[ip[rxn_id, 2]] += 1
+            self.stoi_matrix = stoi_matrix
+
         self.species_on_surface = species_on_surface
         self.reactions_on_surface = reactions_on_surface
 
@@ -378,9 +433,10 @@ cdef class SurfaceReactor(ReactionSystem):
         cdef np.ndarray[np.float64_t, ndim=2] jacobian, dgdk
         cdef list list_of_coverage_deps
         cdef double surface_site_fraction, total_sites, a, m, E
-
-
-
+        cdef np.ndarray[np.float64_t, ndim=1] coverages, coverages_squared, temperature_scaled_coverages
+        cdef np.ndarray[np.float64_t, ndim=2] thermo_dep_coverage
+        cdef np.ndarray[np.float64_t, ndim=1] free_energy_coverage_corrections, rxns_free_energy_change, corrected_K_eq
+        cdef double sp_free_energy_correction
         ir = self.reactant_indices
         ip = self.product_indices
         equilibrium_constants = self.Keq
@@ -414,14 +470,33 @@ cdef class SurfaceReactor(ReactionSystem):
         V = self.V  # constant volume reactor
         A = self.V * surface_volume_ratio_si  # area
         total_sites = self.surface_site_density.value_si * A  # todo: double check units
-
         for j in range(num_core_species):
             if species_on_surface[j]:
                 C[j] = (N[j] / V) / surface_volume_ratio_si
             else:
                 C[j] = N[j] / V
             #: surface species are in mol/m2, gas phase are in mol/m3
             core_species_concentrations[j] = C[j]
+
+        # 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)
+            thermo_dep_coverage[0, :] = coverages
+            thermo_dep_coverage[1, :] = coverages_squared
+            thermo_dep_coverage[2, :] = coverages_squared * coverages
+            thermo_dep_coverage[3, :] = temperature_scaled_coverages
+            thermo_dep_coverage[4, :] = temperature_scaled_coverages * coverages
+            thermo_dep_coverage[5, :] = temperature_scaled_coverages * coverages_squared
+            free_energy_coverage_corrections = np.empty(len(self.thermo_coeff_matrix), dtype=np.float64)
+            for i, matrix in enumerate(self.thermo_coeff_matrix):
+                free_energy_coverage_corrections[i] = np.diag(np.dot(matrix, thermo_dep_coverage)).sum()
+            rxns_free_energy_change = np.matmul(self.stoi_matrix,free_energy_coverage_corrections)
+            corrected_K_eq = copy.deepcopy(self.Keq)
+            corrected_K_eq *= np.exp(-1 * rxns_free_energy_change / (constants.R * self.T.value_si))
+            kr = kf / corrected_K_eq
 
         # Coverage dependence
         coverage_corrections = np.ones_like(kf, float)

rmgpy/thermo/nasa.pxd

@@ -56,6 +56,7 @@ cdef class NASAPolynomial(HeatCapacityModel):
 cdef class NASA(HeatCapacityModel):
 
     cdef public NASAPolynomial poly1, poly2, poly3
+    cdef public dict _thermo_coverage_dependence
 
     cpdef NASAPolynomial select_polynomial(self, double T)