Merge pull request #2646 for coverage-dependent adsorbate thermochemistry

Coverage-dependent thermochemistry for catalysis

This PR adds in coverage dependent thermodynamic models for heterogeneous catalysis modeling. An adsorbate's enthalpy and entropy should be affected by other adsorbates around it. The PR enables RMG to use a polynomial model to estimate the change of an adsorbate's enthalpy or entropy based on the coverage of other adsorbates on the catalyst surface. By doing this, a Cantera yaml file with thermodynamic coverage dependent data can be generated at the end of a RMG simulation, and it can be used by Cantera (>=3.0) to run PFR simulation. RMG looks into the database to read the thermo coverage dependent models.

Author: rwest

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)