From 69f021cf9467180f4dce96fd1d73b5284a3a1ab9 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Sun, 10 Mar 2024 16:49:42 -0400 Subject: [PATCH 01/49] Add thermodynamic coverage dependent models in surface reactor --- rmgpy/solver/surface.pyx | 80 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 79 insertions(+), 1 deletion(-) diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index 105a1d68f7a..3c20de6ab93 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -84,6 +84,7 @@ cdef class SurfaceReactor(ReactionSystem): sensitivity_threshold=1e-3, sens_conditions=None, coverage_dependence=False, + thermo_coverage_dependence=False, ): ReactionSystem.__init__(self, termination, @@ -103,12 +104,14 @@ 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 self.n_sims = n_sims self.coverage_dependencies = {} + self.thermo_coverage_dependencies = {} def convert_initial_keys_to_species_objects(self, species_dict): """ @@ -195,6 +198,31 @@ 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: + # think about the data structure of thermo coverage dependence + # should not be a list, a dictionary would be better + for spec, parameters in sp.thermo.coverage_dependence.items(): + species_index = self.species_index[spec] + try: + list_of_thermo_coverage_deps = self.thermo_coverage_dependencies[species_index] + except KeyError: # doesn't exist yet + list_of_thermo_coverage_deps = [] + self.thermo_coverage_dependencies[sp_index] = list_of_thermo_coverage_deps + # need to specify the entropy and enthalpy models + # linear, piecewise linear, polynomial, interpolative models + list_of_thermo_coverage_deps.append((sp_index, parameters)) + + """ + self.thermo_coverage_dependencies[2] = [(3, {"model":"linear","enthalpy-coefficients":[], "entropy-coefficients":[]}),] + self.thermo_coverage_dependencies[2] = [(3, {"model":"polynomial","enthalpy-coefficients":[], "entropy-coefficients":[]}),] + self.thermo_coverage_dependencies[2] = [(3, {"model":"piecewise-linear","enthalpy-coefficients":{"enthalpy_high":, "enthalpy_low":, "enthalpy_change":}, + "entropy-coefficients":{"entropy_high":, "entropy_low":, "entropy_change"}, "Cp":{"heat-capacity-a":, "heat-capacity-b":}}),] + self.thermo_coverage_dependencies[2] = {(3, {"model":"interpolative","enthalpy-coefficients":{"enthalpy-coverages":, "enthalpies":}, "entropy-coefficients":{"entropy-coverages":, "entropies":}}),} + means that Species with index 2 in the current simulation is used in + Species 3 with parameters for linear, polynomial, piecewise-linear, and interpolative models + """ self.species_on_surface = species_on_surface self.reactions_on_surface = reactions_on_surface @@ -422,7 +450,57 @@ cdef class SurfaceReactor(ReactionSystem): 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 + free_energy_coverage_corrections = np.zeros(len(self.sp_index), float) # length of core + edge species + if self.thermo_coverage_dependence: + """ + self.thermo_coverage_dependencies[2] = [(3, {"model":"linear","enthalpy-coefficients":[], "entropy-coefficients":[]}),] + self.thermo_coverage_dependencies[2] = [(3, {"model":"polynomial","enthalpy-coefficients":[], "entropy-coefficients":[]}),] + self.thermo_coverage_dependencies[2] = [(3, {"model":"piecewise-linear","enthalpy-coefficients":{"enthalpy_high":, "enthalpy_low":, "enthalpy_change":}, + "entropy-coefficients":{"entropy_high":, "entropy_low":, "entropy_change"}, "Cp":{"heat-capacity-a":, "heat-capacity-b":}}),] + self.thermo_coverage_dependencies[2] = {(3, {"model":"interpolative","enthalpy-coefficients":{"enthalpy-coverages":, "enthalpies":}, "entropy-coefficients":{"entropy-coverages":, "entropies":}}),} + means that Species with index 2 in the current simulation is used in + Species 3 with parameters for linear, polynomial, piecewise-linear, and interpolative models + """ + for i, list_of_thermo_coverage_deps in self.thermo_coverage_dependencies.items(): + surface_site_fraction = N[i] / total_sites + if surface_site_fraction < 1e-15: + continue + for j, parameters in list_of_thermo_coverage_deps: + # Species i, Species j + # need to specify the entropy and enthalpy models + # linear, piecewise linear, polynomial, interpolative models + if parameters['model'] == "linear": + pass + elif parameters['model'] == "polynomial": + enthalpy_cov_correction = np.polynomial.polynomial.polyval(surface_site_fraction, parameters['enthalpy-coefficients'].insert(0,0)) # insert 0 for the constant term + entropy_cov_correction = np.polynomial.polynomial.polyval(surface_site_fraction, parameters['entropy-coefficients'].insert(0,0)) + free_energy_coverage_corrections[j] += enthalpy_cov_correction - self.T.value_si * entropy_cov_correction + elif parameters['model'] == "piecewise-linear": + pass + elif parameters['model'] == "interpolative": + pass + corrected_K_eq = copy.deepcopy(self.K_eq) + # correct the K_eq + for j in range(ir.shape[0]): + if ir[j, 0] >= num_core_species or ir[j, 1] >= num_core_species or ir[j, 2] >= num_core_species: + pass + elif ir[j, 1] == -1: # only one reactant + corrected_K_eq[j] *= np.exp(free_energy_coverage_corrections[ir[j, 0]] / (constants.R * self.T.value_si)) + elif ir[j, 2] == -1: # only two reactants + corrected_K_eq[j] *= np.exp((free_energy_coverage_corrections[ir[j, 0]] + free_energy_coverage_corrections[ir[j, 1]]) / (constants.R * self.T.value_si)) + else: # three reactants!! (really?) + corrected_K_eq[j] *= np.exp((free_energy_coverage_corrections[ir[j, 0]] + free_energy_coverage_corrections[ir[j, 1]] + free_energy_coverage_corrections[ir[j, 2]]) / (constants.R * self.T.value_si)) + if ip[j, 0] >= num_core_species or ip[j, 1] >= num_core_species or ip[j, 2] >= num_core_species: + pass + elif ip[j, 1] == -1: # only one reactant + corrected_K_eq[j] /= np.exp(free_energy_coverage_corrections[ir[j, 0]] / (constants.R * self.T.value_si)) + elif ip[j, 2] == -1: # only two reactants + corrected_K_eq[j] /= np.exp((free_energy_coverage_corrections[ir[j, 0]] + free_energy_coverage_corrections[ir[j, 1]]) / (constants.R * self.T.value_si)) + else: # three reactants!! (really?) + corrected_K_eq[j] /= np.exp((free_energy_coverage_corrections[ir[j, 0]] + free_energy_coverage_corrections[ir[j, 1]] + free_energy_coverage_corrections[ir[j, 2]]) / (constants.R * self.T.value_si)) + kr = kf / corrected_K_eq # Coverage dependence coverage_corrections = np.ones_like(kf, float) if self.coverage_dependence: From 1ad7c77ff7602f18385fa37f248cbecb0d536703 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 11 Mar 2024 13:26:58 -0400 Subject: [PATCH 02/49] add thermo coverage dependence instance while reading the file --- rmgpy/rmg/input.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index 04f72c0fa79..f4780e28433 100644 --- a/rmgpy/rmg/input.py +++ b/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)) From 416082c1759d13aa3623eea4f42080cf3a4b9579 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 11 Mar 2024 13:27:21 -0400 Subject: [PATCH 03/49] add import copy package --- rmgpy/solver/surface.pyx | 1 + 1 file changed, 1 insertion(+) diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index 3c20de6ab93..eb9f2b09960 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -43,6 +43,7 @@ cimport rmgpy.constants as constants from rmgpy.quantity import Quantity from rmgpy.quantity cimport ScalarQuantity from rmgpy.solver.base cimport ReactionSystem +import copy cdef class SurfaceReactor(ReactionSystem): """ From 09456561a52a321c8f396527a0995c01bd37c6e2 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 11 Mar 2024 15:30:01 -0400 Subject: [PATCH 04/49] add the instances for thermo coverage dependence --- rmgpy/solver/surface.pyx | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index eb9f2b09960..54151c49938 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -70,6 +70,9 @@ cdef class SurfaceReactor(ReactionSystem): cdef public bint coverage_dependence cdef public dict coverage_dependencies + cdef public bint thermo_coverage_dependence + cdef public dict thermo_coverage_dependencies + def __init__(self, @@ -453,7 +456,7 @@ cdef class SurfaceReactor(ReactionSystem): core_species_concentrations[j] = C[j] # Thermodynamic coverage dependence - free_energy_coverage_corrections = np.zeros(len(self.sp_index), float) # length of core + edge species + free_energy_coverage_corrections = np.zeros(len(self.species_index), float) # length of core + edge species if self.thermo_coverage_dependence: """ self.thermo_coverage_dependencies[2] = [(3, {"model":"linear","enthalpy-coefficients":[], "entropy-coefficients":[]}),] From a42f58fe9114502c9b5a9e26bc539d26cedb0bb8 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 11 Mar 2024 20:43:28 -0400 Subject: [PATCH 05/49] add polynomial thermodynamic coverage model to Wilhoit, NASA, and thermodata --- rmgpy/thermo/nasa.pyx | 42 ++++++++++++++++++++--------- rmgpy/thermo/thermodata.pyx | 52 +++++++++++++++++++++++------------ rmgpy/thermo/wilhoit.pyx | 54 +++++++++++++++++++++++++------------ 3 files changed, 102 insertions(+), 46 deletions(-) diff --git a/rmgpy/thermo/nasa.pyx b/rmgpy/thermo/nasa.pyx index e484f3036f3..832183ac757 100644 --- a/rmgpy/thermo/nasa.pyx +++ b/rmgpy/thermo/nasa.pyx @@ -43,19 +43,20 @@ cdef class NASAPolynomial(HeatCapacityModel): seven-coefficient and nine-coefficient variations are supported. The attributes are: - =============== ============================================================ - Attribute Description - =============== ============================================================ - `coeffs` The seven or nine NASA polynomial coefficients - `Tmin` The minimum temperature in K at which the model is valid, or zero if unknown or undefined - `Tmax` The maximum temperature in K at which the model is valid, or zero if unknown or undefined - `E0` The energy at zero Kelvin (including zero point energy) - `comment` Information about the model (e.g. its source) - =============== ============================================================ + =========================== ============================================================ + Attribute Description + =========================== ============================================================ + `coeffs` The seven or nine NASA polynomial coefficients + `Tmin` The minimum temperature in K at which the model is valid, or zero if unknown or undefined + `Tmax` The maximum temperature in K at which the model is valid, or zero if unknown or undefined + `E0` The energy at zero Kelvin (including zero point energy) + `thermo_coverage_dependence` The coverage dependence of the thermo + `comment` Information about the model (e.g. its source) + =========================== ============================================================ """ - def __init__(self, coeffs=None, Tmin=None, Tmax=None, E0=None, label='', comment=''): + def __init__(self, coeffs=None, Tmin=None, Tmax=None, E0=None, label='', thermo_coverage_dependence=None, comment=''): HeatCapacityModel.__init__(self, Tmin=Tmin, Tmax=Tmax, E0=E0, label=label, comment=comment) self.coeffs = coeffs @@ -73,6 +74,7 @@ cdef class NASAPolynomial(HeatCapacityModel): if self.Tmax is not None: string += ', Tmax={0!r}'.format(self.Tmax) if self.E0 is not None: string += ', E0={0!r}'.format(self.E0) if self.label != '': string += ', label="""{0}"""'.format(self.label) + if self.thermo_coverage_dependence is not None: string += ', thermo_coverage_dependence={0!r}'.format(self.thermo_coverage_dependence) if self.comment != '': string += ', comment="""{0}"""'.format(self.comment) string += ')' return string @@ -81,7 +83,7 @@ cdef class NASAPolynomial(HeatCapacityModel): """ A helper function used when pickling an object. """ - return (NASAPolynomial, ([self.cm2, self.cm1, self.c0, self.c1, self.c2, self.c3, self.c4, self.c5, self.c6], self.Tmin, self.Tmax, self.E0, self.label, self.comment)) + return (NASAPolynomial, ([self.cm2, self.cm1, self.c0, self.c1, self.c2, self.c3, self.c4, self.c5, self.c6], self.Tmin, self.Tmax, self.E0, self.label, self.thermo_coverage_dependence, self.comment)) cpdef dict as_dict(self): output_dictionary = super(NASAPolynomial, self).as_dict() @@ -108,6 +110,21 @@ cdef class NASAPolynomial(HeatCapacityModel): else: raise ValueError('Invalid number of NASA polynomial coefficients; expected 7 or 9, got {0:d}.'.format(len(value))) + property thermo_coverage_dependence: + """The coverage dependence of the thermo""" + def __get__(self): + return self._thermo_coverage_dependence + def __set__(self, value): + self._thermo_coverage_dependence = {} + if value: + for species, parameters in value.items(): + # just the polynomial model for now + processed_parameters = {'model': parameters['model'], + 'enthalpy-coefficients': [quantity.Dimensionless(p) for p in parameters['enthalpy-coefficients']], + 'entropy-coefficients': [quantity.Dimensionless(p) for p in parameters['entropy-coefficients']], + } + self._coverage_dependence[species] = processed_parameters + cpdef double get_heat_capacity(self, double T) except -1000000000: """ Return the constant-pressure heat capacity in J/mol*K at the specified @@ -347,6 +364,7 @@ cdef class NASA(HeatCapacityModel): Cp0 = self.Cp0, CpInf = self.CpInf, E0 = self.E0, + thermo_coverage_dependence = self.thermo_coverage_dependence, comment = self.comment ) @@ -379,7 +397,7 @@ cdef class NASA(HeatCapacityModel): H298 = self.get_enthalpy(298) S298 = self.get_entropy(298) - return Wilhoit(label=self.label,comment=self.comment).fit_to_data(Tdata, Cpdata, Cp0, CpInf, H298, S298) + return Wilhoit(label=self.label, thermo_coverage_dependence=self.thermo_coverage_dependence, comment=self.comment).fit_to_data(Tdata, Cpdata, Cp0, CpInf, H298, S298) cpdef NASA change_base_enthalpy(self, double deltaH): """ diff --git a/rmgpy/thermo/thermodata.pyx b/rmgpy/thermo/thermodata.pyx index 5248e79a563..99c72377897 100644 --- a/rmgpy/thermo/thermodata.pyx +++ b/rmgpy/thermo/thermodata.pyx @@ -43,27 +43,29 @@ cdef class ThermoData(HeatCapacityModel): A heat capacity model based on a set of discrete heat capacity data points. The attributes are: - =============== ============================================================ - Attribute Description - =============== ============================================================ - `Tdata` An array of temperatures at which the heat capacity is known - `Cpdata` An array of heat capacities at the given temperatures - `H298` The standard enthalpy of formation at 298 K - `S298` The standard entropy at 298 K - `Tmin` The minimum temperature at which the model is valid, or zero if unknown or undefined - `Tmax` The maximum temperature at which the model is valid, or zero if unknown or undefined - `E0` The energy at zero Kelvin (including zero point energy) - `comment` Information about the model (e.g. its source) - =============== ============================================================ + =========================== ============================================================ + Attribute Description + =========================== ============================================================ + `Tdata` An array of temperatures at which the heat capacity is known + `Cpdata` An array of heat capacities at the given temperatures + `H298` The standard enthalpy of formation at 298 K + `S298` The standard entropy at 298 K + `Tmin` The minimum temperature at which the model is valid, or zero if unknown or undefined + `Tmax` The maximum temperature at which the model is valid, or zero if unknown or undefined + `E0` The energy at zero Kelvin (including zero point energy) + `thermo_coverage_dependence` The coverage dependence of the thermo + `comment` Information about the model (e.g. its source) + =========================== ============================================================ """ - def __init__(self, Tdata=None, Cpdata=None, H298=None, S298=None, Cp0=None, CpInf=None, Tmin=None, Tmax=None, E0=None, label = '',comment=''): + def __init__(self, Tdata=None, Cpdata=None, H298=None, S298=None, Cp0=None, CpInf=None, Tmin=None, Tmax=None, E0=None, label = '', thermo_coverage_dependence=None, comment=''): HeatCapacityModel.__init__(self, Tmin=Tmin, Tmax=Tmax, E0=E0, Cp0=Cp0, CpInf=CpInf, label=label, comment=comment) self.H298 = H298 self.S298 = S298 self.Tdata = Tdata self.Cpdata = Cpdata + self.thermo_coverage_dependence = thermo_coverage_dependence def __repr__(self): """ @@ -77,6 +79,7 @@ cdef class ThermoData(HeatCapacityModel): if self.Tmax is not None: string += ', Tmax={0!r}'.format(self.Tmax) if self.E0 is not None: string += ', E0={0!r}'.format(self.E0) if self.label != '': string +=', label="""{0}"""'.format(self.label) + if self.thermo_coverage_dependence is not None: string += ', thermo_coverage_dependence={0!r}'.format(self.thermo_coverage_dependence) if self.comment != '': string += ', comment="""{0}"""'.format(self.comment) string += ')' return string @@ -85,7 +88,7 @@ cdef class ThermoData(HeatCapacityModel): """ A helper function used when pickling a ThermoData object. """ - return (ThermoData, (self.Tdata, self.Cpdata, self.H298, self.S298, self.Cp0, self.CpInf, self.Tmin, self.Tmax, self.E0, self.label, self.comment)) + return (ThermoData, (self.Tdata, self.Cpdata, self.H298, self.S298, self.Cp0, self.CpInf, self.Tmin, self.Tmax, self.E0, self.label, self.thermo_coverage_dependence, self.comment)) property Tdata: """An array of temperatures at which the heat capacity is known.""" @@ -114,6 +117,21 @@ cdef class ThermoData(HeatCapacityModel): return self._S298 def __set__(self, value): self._S298 = quantity.Entropy(value) + + property thermo_coverage_dependence: + """The coverage dependence of the thermo""" + def __get__(self): + return self._thermo_coverage_dependence + def __set__(self, value): + self._thermo_coverage_dependence = {} + if value: + for species, parameters in value.items(): + # just the polynomial model for now + processed_parameters = {'model': parameters['model'], + 'enthalpy-coefficients': [quantity.Dimensionless(p) for p in parameters['enthalpy-coefficients']], + 'entropy-coefficients': [quantity.Dimensionless(p) for p in parameters['entropy-coefficients']], + } + self._coverage_dependence[species] = processed_parameters @cython.boundscheck(False) @cython.wraparound(False) @@ -361,12 +379,12 @@ cdef class ThermoData(HeatCapacityModel): S298 = self.get_entropy(298) Cp0 = self._Cp0.value_si CpInf = self._CpInf.value_si - + thermo_cov_dep = self.thermo_coverage_dependence if B: - return Wilhoit(label=self.label,comment=self.comment).fit_to_data_for_constant_b(Tdata, Cpdata, Cp0, CpInf, H298, S298, B=B) + return Wilhoit(label=self.label, thermo_coverage_dependence=thermo_cov_dep, comment=self.comment).fit_to_data_for_constant_b(Tdata, Cpdata, Cp0, CpInf, H298, S298, B=B) else: - return Wilhoit(label=self.label,comment=self.comment).fit_to_data(Tdata, Cpdata, Cp0, CpInf, H298, S298) + return Wilhoit(label=self.label, thermo_coverage_dependence=thermo_cov_dep, comment=self.comment).fit_to_data(Tdata, Cpdata, Cp0, CpInf, H298, S298) cpdef NASA to_nasa(self, double Tmin, double Tmax, double Tint, bint fixedTint=False, bint weighting=True, int continuity=3): """ diff --git a/rmgpy/thermo/wilhoit.pyx b/rmgpy/thermo/wilhoit.pyx index 20a1061b18e..2dc974b016d 100644 --- a/rmgpy/thermo/wilhoit.pyx +++ b/rmgpy/thermo/wilhoit.pyx @@ -45,24 +45,25 @@ cdef class Wilhoit(HeatCapacityModel): """ A heat capacity model based on the Wilhoit equation. The attributes are: - =============== ============================================================ - Attribute Description - =============== ============================================================ - `a0` The zeroth-order Wilhoit polynomial coefficient - `a1` The first-order Wilhoit polynomial coefficient - `a2` The second-order Wilhoit polynomial coefficient - `a3` The third-order Wilhoit polynomial coefficient - `H0` The integration constant for enthalpy (not H at T=0) - `S0` The integration constant for entropy (not S at T=0) - `E0` The energy at zero Kelvin (including zero point energy) - `B` The Wilhoit scaled temperature coefficient in K - `Tmin` The minimum temperature in K at which the model is valid, or zero if unknown or undefined - `Tmax` The maximum temperature in K at which the model is valid, or zero if unknown or undefined - `comment` Information about the model (e.g. its source) - =============== ============================================================ + ============================= ========================================================================================= + Attribute Description + ============================= ========================================================================================= + `a0` The zeroth-order Wilhoit polynomial coefficient + `a1` The first-order Wilhoit polynomial coefficient + `a2` The second-order Wilhoit polynomial coefficient + `a3` The third-order Wilhoit polynomial coefficient + `H0` The integration constant for enthalpy (not H at T=0) + `S0` The integration constant for entropy (not S at T=0) + `E0` The energy at zero Kelvin (including zero point energy) + `B` The Wilhoit scaled temperature coefficient in K + `Tmin` The minimum temperature in K at which the model is valid, or zero if unknown or undefined + `Tmax` The maximum temperature in K at which the model is valid, or zero if unknown or undefined + `thermo_covreage_dependence` The coverage dependence of the thermo + `comment` Information about the model (e.g. its source) + ============================= ========================================================================================= """ - def __init__(self, Cp0=None, CpInf=None, a0=0.0, a1=0.0, a2=0.0, a3=0.0, H0=None, S0=None, B=None, Tmin=None, Tmax=None, label='', comment=''): + def __init__(self, Cp0=None, CpInf=None, a0=0.0, a1=0.0, a2=0.0, a3=0.0, H0=None, S0=None, B=None, Tmin=None, Tmax=None, label='', thermo_coverage_dependence=None, comment=''): HeatCapacityModel.__init__(self, Tmin=Tmin, Tmax=Tmax, Cp0=Cp0, CpInf=CpInf, label=label, comment=comment) self.B = B self.a0 = a0 @@ -71,6 +72,7 @@ cdef class Wilhoit(HeatCapacityModel): self.a3 = a3 self.H0 = H0 self.S0 = S0 + self.thermo_coverage_dependence = thermo_coverage_dependence def __repr__(self): """ @@ -82,6 +84,7 @@ cdef class Wilhoit(HeatCapacityModel): if self.Tmin is not None: string += ', Tmin={0!r}'.format(self.Tmin) if self.Tmax is not None: string += ', Tmax={0!r}'.format(self.Tmax) if self.label != '': string += ', label="""{0}"""'.format(self.label) + if self.thermo_coverage_dependence is not None: string += ', thermo_coverage_dependence={0!r}'.format(self.thermo_coverage_dependence) if self.comment != '': string += ', comment="""{0}"""'.format(self.comment) string += ')' return string @@ -90,7 +93,7 @@ cdef class Wilhoit(HeatCapacityModel): """ A helper function used when pickling a Wilhoit object. """ - return (Wilhoit, (self.Cp0, self.CpInf, self.a0, self.a1, self.a2, self.a3, self.H0, self.S0, self.B, self.Tmin, self.Tmax, self.label, self.comment)) + return (Wilhoit, (self.Cp0, self.CpInf, self.a0, self.a1, self.a2, self.a3, self.H0, self.S0, self.B, self.Tmin, self.Tmax, self.label, self.thermo_coverage_dependence, self.comment)) cpdef dict as_dict(self): """ @@ -136,6 +139,21 @@ cdef class Wilhoit(HeatCapacityModel): return self._S0 def __set__(self, value): self._S0 = quantity.Entropy(value) + + property thermo_coverage_dependence: + """The coverage dependence of the thermo""" + def __get__(self): + return self._thermo_coverage_dependence + def __set__(self, value): + self._thermo_coverage_dependence = {} + if value: + for species, parameters in value.items(): + # just the polynomial model for now + processed_parameters = {'model': parameters['model'], + 'enthalpy-coefficients': [quantity.Dimensionless(p) for p in parameters['enthalpy-coefficients']], + 'entropy-coefficients': [quantity.Dimensionless(p) for p in parameters['entropy-coefficients']], + } + self._coverage_dependence[species] = processed_parameters cpdef double get_heat_capacity(self, double T) except -1000000000: """ @@ -478,6 +496,7 @@ cdef class Wilhoit(HeatCapacityModel): Cp0 = self.Cp0, CpInf = self.CpInf, E0 = self.E0, + thermo_coverage_dependence = self.thermo_coverage_dependence, comment = self.comment ) @@ -582,6 +601,7 @@ cdef class Wilhoit(HeatCapacityModel): Cp0 = self.Cp0, CpInf = self.CpInf, label = self.label, + thermo_coverage_dependence = self.thermo_coverage_dependence, comment = self.comment, ) From 6d4ab278989c020e249fdc9ed40bd52c6f3c1630 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Tue, 12 Mar 2024 14:47:25 -0400 Subject: [PATCH 06/49] declare _thermo_coverage_dependence in the pxd files --- rmgpy/thermo/nasa.pxd | 1 + rmgpy/thermo/thermodata.pxd | 1 + rmgpy/thermo/wilhoit.pxd | 1 + 3 files changed, 3 insertions(+) diff --git a/rmgpy/thermo/nasa.pxd b/rmgpy/thermo/nasa.pxd index b462301f72a..fb221419dae 100644 --- a/rmgpy/thermo/nasa.pxd +++ b/rmgpy/thermo/nasa.pxd @@ -34,6 +34,7 @@ from rmgpy.thermo.wilhoit cimport Wilhoit cdef class NASAPolynomial(HeatCapacityModel): cdef public double cm2, cm1, c0, c1, c2, c3, c4, c5, c6 + cdef public dict _thermo_coverage_dependence cpdef dict as_dict(self) diff --git a/rmgpy/thermo/thermodata.pxd b/rmgpy/thermo/thermodata.pxd index d9052db0b4e..0a1d08b12ff 100644 --- a/rmgpy/thermo/thermodata.pxd +++ b/rmgpy/thermo/thermodata.pxd @@ -38,6 +38,7 @@ cdef class ThermoData(HeatCapacityModel): cdef public ScalarQuantity _H298, _S298 cdef public ArrayQuantity _Tdata, _Cpdata + cdef public dict _thermo_coverage_dependence cpdef double get_heat_capacity(self, double T) except -1000000000 diff --git a/rmgpy/thermo/wilhoit.pxd b/rmgpy/thermo/wilhoit.pxd index 576c8a2e671..e2cd2419e42 100644 --- a/rmgpy/thermo/wilhoit.pxd +++ b/rmgpy/thermo/wilhoit.pxd @@ -38,6 +38,7 @@ cdef class Wilhoit(HeatCapacityModel): cdef public ScalarQuantity _B, _H0, _S0 cdef public double a0, a1, a2, a3 + cdef public dict _thermo_coverage_dependence cpdef dict as_dict(self) From f23a8d600fe6e5ea58f9eb4d8bece73d5d13701b Mon Sep 17 00:00:00 2001 From: 12Chao Date: Tue, 12 Mar 2024 14:48:37 -0400 Subject: [PATCH 07/49] fix the typos, and add the missing import command line --- rmgpy/thermo/nasa.pyx | 3 ++- rmgpy/thermo/thermodata.pyx | 2 +- rmgpy/thermo/wilhoit.pyx | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/rmgpy/thermo/nasa.pyx b/rmgpy/thermo/nasa.pyx index 832183ac757..8b2a65ae946 100644 --- a/rmgpy/thermo/nasa.pyx +++ b/rmgpy/thermo/nasa.pyx @@ -34,6 +34,7 @@ cimport numpy as np from libc.math cimport log cimport rmgpy.constants as constants +import rmgpy.quantity as quantity ################################################################################ @@ -123,7 +124,7 @@ cdef class NASAPolynomial(HeatCapacityModel): 'enthalpy-coefficients': [quantity.Dimensionless(p) for p in parameters['enthalpy-coefficients']], 'entropy-coefficients': [quantity.Dimensionless(p) for p in parameters['entropy-coefficients']], } - self._coverage_dependence[species] = processed_parameters + self._thermo_coverage_dependence[species] = processed_parameters cpdef double get_heat_capacity(self, double T) except -1000000000: """ diff --git a/rmgpy/thermo/thermodata.pyx b/rmgpy/thermo/thermodata.pyx index 99c72377897..e38ae87b51b 100644 --- a/rmgpy/thermo/thermodata.pyx +++ b/rmgpy/thermo/thermodata.pyx @@ -131,7 +131,7 @@ cdef class ThermoData(HeatCapacityModel): 'enthalpy-coefficients': [quantity.Dimensionless(p) for p in parameters['enthalpy-coefficients']], 'entropy-coefficients': [quantity.Dimensionless(p) for p in parameters['entropy-coefficients']], } - self._coverage_dependence[species] = processed_parameters + self._thermo_coverage_dependence[species] = processed_parameters @cython.boundscheck(False) @cython.wraparound(False) diff --git a/rmgpy/thermo/wilhoit.pyx b/rmgpy/thermo/wilhoit.pyx index 2dc974b016d..15299666364 100644 --- a/rmgpy/thermo/wilhoit.pyx +++ b/rmgpy/thermo/wilhoit.pyx @@ -153,7 +153,7 @@ cdef class Wilhoit(HeatCapacityModel): 'enthalpy-coefficients': [quantity.Dimensionless(p) for p in parameters['enthalpy-coefficients']], 'entropy-coefficients': [quantity.Dimensionless(p) for p in parameters['entropy-coefficients']], } - self._coverage_dependence[species] = processed_parameters + self._thermo_coverage_dependence[species] = processed_parameters cpdef double get_heat_capacity(self, double T) except -1000000000: """ From 50402b8e0441913b0c3aaf0594fe8606135fba8a Mon Sep 17 00:00:00 2001 From: 12Chao Date: Tue, 12 Mar 2024 16:12:12 -0400 Subject: [PATCH 08/49] Add thermo_coverage_dependence as new property to NASA --- rmgpy/thermo/nasa.pxd | 1 + rmgpy/thermo/nasa.pyx | 36 ++++++++++++++++++++++++++---------- 2 files changed, 27 insertions(+), 10 deletions(-) diff --git a/rmgpy/thermo/nasa.pxd b/rmgpy/thermo/nasa.pxd index fb221419dae..43f629e60c9 100644 --- a/rmgpy/thermo/nasa.pxd +++ b/rmgpy/thermo/nasa.pxd @@ -57,6 +57,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) diff --git a/rmgpy/thermo/nasa.pyx b/rmgpy/thermo/nasa.pyx index 8b2a65ae946..d7af53acd06 100644 --- a/rmgpy/thermo/nasa.pyx +++ b/rmgpy/thermo/nasa.pyx @@ -221,19 +221,20 @@ cdef class NASA(HeatCapacityModel): A heat capacity model based on a set of one, two, or three :class:`NASAPolynomial` objects. The attributes are: - =============== ============================================================ - Attribute Description - =============== ============================================================ - `polynomials` The list of NASA polynomials to use in this model - `Tmin` The minimum temperature in K at which the model is valid, or zero if unknown or undefined - `Tmax` The maximum temperature in K at which the model is valid, or zero if unknown or undefined - `E0` The energy at zero Kelvin (including zero point energy) - `comment` Information about the model (e.g. its source) - =============== ============================================================ + ============================ ============================================================ + Attribute Description + ============================ ============================================================ + `polynomials` The list of NASA polynomials to use in this model + `Tmin` The minimum temperature in K at which the model is valid, or zero if unknown or undefined + `Tmax` The maximum temperature in K at which the model is valid, or zero if unknown or undefined + `E0` The energy at zero Kelvin (including zero point energy) + `thermo_coverage_dependence` The coverage dependence of the thermo + `comment` Information about the model (e.g. its source) + ============================ ============================================================ """ - def __init__(self, polynomials=None, Tmin=None, Tmax=None, E0=None, Cp0=None, CpInf=None, label='', comment=''): + def __init__(self, polynomials=None, Tmin=None, Tmax=None, E0=None, Cp0=None, CpInf=None, label='', thermo_coverage_dependence=None, comment=''): HeatCapacityModel.__init__(self, Tmin=Tmin, Tmax=Tmax, E0=E0, Cp0=Cp0, CpInf=CpInf, label=label, comment=comment) if isinstance(polynomials, dict): self.polynomials = list(polynomials.values()) @@ -318,6 +319,21 @@ cdef class NASA(HeatCapacityModel): else: raise ValueError('No valid NASA polynomial at temperature {0:g} K.'.format(T)) + property thermo_coverage_dependence: + """The coverage dependence of the thermo""" + def __get__(self): + return self._thermo_coverage_dependence + def __set__(self, value): + self._thermo_coverage_dependence = {} + if value: + for species, parameters in value.items(): + # just the polynomial model for now + processed_parameters = {'model': parameters['model'], + 'enthalpy-coefficients': [quantity.Dimensionless(p) for p in parameters['enthalpy-coefficients']], + 'entropy-coefficients': [quantity.Dimensionless(p) for p in parameters['entropy-coefficients']], + } + self._thermo_coverage_dependence[species] = processed_parameters + cpdef double get_heat_capacity(self, double T) except -1000000000: """ Return the constant-pressure heat capacity From 669ba4ee39e6711c5ac61dd388a6accf03ae0ca4 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Wed, 13 Mar 2024 14:53:34 -0400 Subject: [PATCH 09/49] Move thermo_coverage_dependence from NASAPolynomial to NASA --- rmgpy/thermo/nasa.pxd | 1 - rmgpy/thermo/nasa.pyx | 44 +++++++++++++++---------------------------- 2 files changed, 15 insertions(+), 30 deletions(-) diff --git a/rmgpy/thermo/nasa.pxd b/rmgpy/thermo/nasa.pxd index 43f629e60c9..6922c51e540 100644 --- a/rmgpy/thermo/nasa.pxd +++ b/rmgpy/thermo/nasa.pxd @@ -34,7 +34,6 @@ from rmgpy.thermo.wilhoit cimport Wilhoit cdef class NASAPolynomial(HeatCapacityModel): cdef public double cm2, cm1, c0, c1, c2, c3, c4, c5, c6 - cdef public dict _thermo_coverage_dependence cpdef dict as_dict(self) diff --git a/rmgpy/thermo/nasa.pyx b/rmgpy/thermo/nasa.pyx index d7af53acd06..4014cdfc1be 100644 --- a/rmgpy/thermo/nasa.pyx +++ b/rmgpy/thermo/nasa.pyx @@ -44,20 +44,19 @@ cdef class NASAPolynomial(HeatCapacityModel): seven-coefficient and nine-coefficient variations are supported. The attributes are: - =========================== ============================================================ - Attribute Description - =========================== ============================================================ - `coeffs` The seven or nine NASA polynomial coefficients - `Tmin` The minimum temperature in K at which the model is valid, or zero if unknown or undefined - `Tmax` The maximum temperature in K at which the model is valid, or zero if unknown or undefined - `E0` The energy at zero Kelvin (including zero point energy) - `thermo_coverage_dependence` The coverage dependence of the thermo - `comment` Information about the model (e.g. its source) - =========================== ============================================================ + ========= ============================================================ + Attribute Description + ========= ============================================================ + `coeffs` The seven or nine NASA polynomial coefficients + `Tmin` The minimum temperature in K at which the model is valid, or zero if unknown or undefined + `Tmax` The maximum temperature in K at which the model is valid, or zero if unknown or undefined + `E0` The energy at zero Kelvin (including zero point energy) + `comment` Information about the model (e.g. its source) + ========= ============================================================ """ - def __init__(self, coeffs=None, Tmin=None, Tmax=None, E0=None, label='', thermo_coverage_dependence=None, comment=''): + def __init__(self, coeffs=None, Tmin=None, Tmax=None, E0=None, label='', comment=''): HeatCapacityModel.__init__(self, Tmin=Tmin, Tmax=Tmax, E0=E0, label=label, comment=comment) self.coeffs = coeffs @@ -75,7 +74,6 @@ cdef class NASAPolynomial(HeatCapacityModel): if self.Tmax is not None: string += ', Tmax={0!r}'.format(self.Tmax) if self.E0 is not None: string += ', E0={0!r}'.format(self.E0) if self.label != '': string += ', label="""{0}"""'.format(self.label) - if self.thermo_coverage_dependence is not None: string += ', thermo_coverage_dependence={0!r}'.format(self.thermo_coverage_dependence) if self.comment != '': string += ', comment="""{0}"""'.format(self.comment) string += ')' return string @@ -84,7 +82,7 @@ cdef class NASAPolynomial(HeatCapacityModel): """ A helper function used when pickling an object. """ - return (NASAPolynomial, ([self.cm2, self.cm1, self.c0, self.c1, self.c2, self.c3, self.c4, self.c5, self.c6], self.Tmin, self.Tmax, self.E0, self.label, self.thermo_coverage_dependence, self.comment)) + return (NASAPolynomial, ([self.cm2, self.cm1, self.c0, self.c1, self.c2, self.c3, self.c4, self.c5, self.c6], self.Tmin, self.Tmax, self.E0, self.label, self.comment)) cpdef dict as_dict(self): output_dictionary = super(NASAPolynomial, self).as_dict() @@ -111,21 +109,6 @@ cdef class NASAPolynomial(HeatCapacityModel): else: raise ValueError('Invalid number of NASA polynomial coefficients; expected 7 or 9, got {0:d}.'.format(len(value))) - property thermo_coverage_dependence: - """The coverage dependence of the thermo""" - def __get__(self): - return self._thermo_coverage_dependence - def __set__(self, value): - self._thermo_coverage_dependence = {} - if value: - for species, parameters in value.items(): - # just the polynomial model for now - processed_parameters = {'model': parameters['model'], - 'enthalpy-coefficients': [quantity.Dimensionless(p) for p in parameters['enthalpy-coefficients']], - 'entropy-coefficients': [quantity.Dimensionless(p) for p in parameters['entropy-coefficients']], - } - self._thermo_coverage_dependence[species] = processed_parameters - cpdef double get_heat_capacity(self, double T) except -1000000000: """ Return the constant-pressure heat capacity in J/mol*K at the specified @@ -240,6 +223,7 @@ cdef class NASA(HeatCapacityModel): self.polynomials = list(polynomials.values()) else: self.polynomials = polynomials + self.thermo_coverage_dependence = thermo_coverage_dependence def __repr__(self): """ @@ -254,6 +238,7 @@ cdef class NASA(HeatCapacityModel): if self.Cp0 is not None: string += ', Cp0={0!r}'.format(self.Cp0) if self.CpInf is not None: string += ', CpInf={0!r}'.format(self.CpInf) if self.label != '': string += ', label="""{0}"""'.format(self.label) + if self.thermo_coverage_dependence is not None: string += ', thermo_coverage_dependence={0!r}'.format(self.thermo_coverage_dependence) if self.comment != '': string += ', comment="""{0}"""'.format(self.comment) string += ')' return string @@ -262,7 +247,7 @@ cdef class NASA(HeatCapacityModel): """ A helper function used when pickling an object. """ - return (NASA, (self.polynomials, self.Tmin, self.Tmax, self.E0, self.Cp0, self.CpInf, self.label, self.comment)) + return (NASA, (self.polynomials, self.Tmin, self.Tmax, self.E0, self.Cp0, self.CpInf, self.label, self.thermo_coverage_dependence, self.comment)) cpdef dict as_dict(self): """ @@ -434,6 +419,7 @@ 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. From f887b57517d479ae2f342f078f8629d85a6cd66c Mon Sep 17 00:00:00 2001 From: 12Chao Date: Wed, 13 Mar 2024 16:01:00 -0400 Subject: [PATCH 10/49] Add unit test for coverage dependent thermodynamics in RMG thermo classes --- test/rmgpy/thermo/convertTest.py | 11 +++++++++++ test/rmgpy/thermo/nasaTest.py | 14 ++++++++++++++ test/rmgpy/thermo/thermodataTest.py | 10 ++++++++++ test/rmgpy/thermo/wilhoitTest.py | 22 +++++++++++++++++++++- 4 files changed, 56 insertions(+), 1 deletion(-) diff --git a/test/rmgpy/thermo/convertTest.py b/test/rmgpy/thermo/convertTest.py index 94637ce8f63..d35f3de4c7f 100644 --- a/test/rmgpy/thermo/convertTest.py +++ b/test/rmgpy/thermo/convertTest.py @@ -57,6 +57,7 @@ def setup_class(self): S0=(-118.46 * constants.R, "J/(mol*K)"), Tmin=(10, "K"), Tmax=(3000, "K"), + thermo_coverage_dependence = {'OX':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}}, comment="C2H6", ) self.nasa = NASA( @@ -93,6 +94,7 @@ def setup_class(self): E0=(-93.6077, "kJ/mol"), Cp0=(4.0 * constants.R, "J/(mol*K)"), CpInf=(21.5 * constants.R, "J/(mol*K)"), + thermo_coverage_dependence = {'OX':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}}, comment="C2H6", ) self.thermodata = ThermoData( @@ -108,6 +110,7 @@ def setup_class(self): Tmin=(10, "K"), Tmax=(3000, "K"), E0=(-93.6077, "kJ/mol"), + thermo_coverage_dependence = {'OX':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}}, comment="C2H6", ) @@ -129,6 +132,7 @@ def test_convert_wilhoit_to_nasa(self): s_nasa = nasa.get_entropy(T) assert abs(s_nasa - s_wilhoit) < 1e0 assert abs(wilhoit.E0.value_si - nasa.E0.value_si) < 1e1 + assert repr(nasa.thermo_coverage_dependence) == repr(wilhoit.thermo_coverage_dependence) def test_convert_wilhoit_to_thermo_data(self): """ @@ -149,6 +153,7 @@ def test_convert_wilhoit_to_thermo_data(self): s_thermodata = thermodata.get_entropy(T) assert round(abs(s_thermodata - s_wilhoit), 4) == 0 assert abs(wilhoit.E0.value_si - thermodata.E0.value_si) < 1e1 + assert repr(thermodata.thermo_coverage_dependence) == repr(wilhoit.thermo_coverage_dependence) def test_convert_nasa_to_wilhoit(self): """ @@ -168,6 +173,7 @@ def test_convert_nasa_to_wilhoit(self): s_nasa = nasa.get_entropy(T) assert abs(s_nasa - s_wilhoit) < 1e0 assert abs(nasa.E0.value_si - wilhoit.E0.value_si) < 2e1 + assert repr(wilhoit.thermo_coverage_dependence) == repr(nasa.thermo_coverage_dependence) def test_convert_nasa_to_thermo_data(self): """ @@ -188,6 +194,7 @@ def test_convert_nasa_to_thermo_data(self): s_nasa = nasa.get_entropy(T) assert round(abs(s_nasa - s_thermodata), 4) == 0 assert abs(nasa.E0.value_si - thermodata.E0.value_si) < 1e1 + assert repr(thermodata.thermo_coverage_dependence) == repr(nasa.thermo_coverage_dependence) def test_convert_thermo_data_to_wilhoit(self): """ @@ -208,6 +215,7 @@ def test_convert_thermo_data_to_wilhoit(self): s_thermodata = thermodata.get_entropy(T) assert round(abs(s_thermodata - s_wilhoit), 3) == 0 assert abs(thermodata.E0.value_si - wilhoit.E0.value_si) < 1e1 + assert repr(wilhoit.thermo_coverage_dependence) == repr(thermodata.thermo_coverage_dependence) def test_convert_thermo_data_to_nasa(self): """ @@ -228,6 +236,7 @@ def test_convert_thermo_data_to_nasa(self): s_nasa = nasa.get_entropy(T) assert abs(s_nasa - s_thermodata) < 1e0 assert abs(thermodata.E0.value_si - nasa.E0.value_si) < 1e1 + assert repr(nasa.thermo_coverage_dependence) == repr(thermodata.thermo_coverage_dependence) def test_wilhoit_nasa_wilhoit(self): """ @@ -248,6 +257,7 @@ def test_wilhoit_nasa_wilhoit(self): s_2 = wilhoit2.get_entropy(T) assert abs(s_1 - s_2) < 1e0 assert abs(wilhoit1.E0.value_si - wilhoit2.E0.value_si) < 1e1 + assert repr(wilhoit1.thermo_coverage_dependence) == repr(wilhoit2.thermo_coverage_dependence) def test_wilhoit_thermo_data_wilhoit(self): """ @@ -268,3 +278,4 @@ def test_wilhoit_thermo_data_wilhoit(self): s_2 = wilhoit2.get_entropy(T) assert abs(s_1 - s_2) < 1e0 assert abs(wilhoit1.E0.value_si - wilhoit2.E0.value_si) < 1e1 + assert repr(wilhoit1.thermo_coverage_dependence) == repr(wilhoit2.thermo_coverage_dependence) diff --git a/test/rmgpy/thermo/nasaTest.py b/test/rmgpy/thermo/nasaTest.py index f08b28e0da1..8945bef5f42 100644 --- a/test/rmgpy/thermo/nasaTest.py +++ b/test/rmgpy/thermo/nasaTest.py @@ -39,6 +39,7 @@ import rmgpy.constants as constants from rmgpy.quantity import ScalarQuantity from rmgpy.thermo.nasa import NASA, NASAPolynomial +from rmgpy.quantity import Dimensionless class TestNASA: @@ -69,6 +70,7 @@ def setup_class(self): self.Tmax = 3000.0 self.Tint = 650.73 self.E0 = -782292.0 # J/mol. + self.thermo_coverage_dependence = {'OX':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}} self.comment = "C2H6" self.nasa = NASA( polynomials=[ @@ -82,6 +84,7 @@ def setup_class(self): Tmin=(self.Tmin, "K"), Tmax=(self.Tmax, "K"), E0=(self.E0, "J/mol"), + thermo_coverage_dependence=self.thermo_coverage_dependence, comment=self.comment, ) @@ -136,6 +139,12 @@ def test_comment(self): Test that the NASA comment property was properly set. """ assert self.nasa.comment == self.comment + + def test_thermo_coverage_dependence(self): + """ + Test that the thermo_coverage_dependence property was properly set. + """ + assert repr(self.nasa.thermo_coverage_dependence) == repr(self.thermo_coverage_dependence) def test_is_temperature_valid(self): """ @@ -263,6 +272,7 @@ def test_pickle(self): assert self.nasa.Tmax.units == nasa.Tmax.units assert self.nasa.E0.value == nasa.E0.value assert self.nasa.E0.units == nasa.E0.units + assert repr(self.nasa.thermo_coverage_dependence) == repr(nasa.thermo_coverage_dependence) assert self.nasa.comment == nasa.comment def test_repr(self): @@ -296,6 +306,7 @@ def test_repr(self): assert self.nasa.Tmax.units == nasa.Tmax.units assert self.nasa.E0.value == nasa.E0.value assert self.nasa.E0.units == nasa.E0.units + assert repr(self.nasa.thermo_coverage_dependence) == repr(nasa.thermo_coverage_dependence) assert self.nasa.comment == nasa.comment def test_to_cantera(self): @@ -326,6 +337,7 @@ def test_to_nasa(self): spc = Species().from_smiles("CC") spc.get_thermo_data() + spc.thermo.thermo_coverage_dependence = self.thermo_coverage_dependence T = 1350.0 # not 298K! @@ -353,6 +365,7 @@ def test_to_nasa(self): assert round(abs(s_nasa - s_td), -1) == 0 assert wilhoit.comment == nasa.comment + assert repr(wilhoit.thermo_coverage_dependence) == repr(nasa.thermo_coverage_dependence) # nasa to wilhoi performed in wilhoitTest @@ -364,6 +377,7 @@ def test_nasa_as_dict_full(self): assert nasa_dict["E0"]["value"] == self.E0 assert nasa_dict["Tmin"]["value"] == self.Tmin assert nasa_dict["Tmax"]["value"] == self.Tmax + assert repr(nasa_dict["thermo_coverage_dependence"]) == "{'OX': {'model': 'polynomial', 'enthalpy-coefficients': [{'class': 'ScalarQuantity', 'value': 1.0}, {'class': 'ScalarQuantity', 'value': 2.0}, {'class': 'ScalarQuantity', 'value': 3.0}], 'entropy-coefficients': [{'class': 'ScalarQuantity', 'value': 1.0}, {'class': 'ScalarQuantity', 'value': 2.0}, {'class': 'ScalarQuantity', 'value': 3.0}]}}" assert nasa_dict["comment"] == self.comment assert tuple(nasa_dict["polynomials"]["polynomial1"]["coeffs"]["object"]) == tuple(self.coeffs_low) assert tuple(nasa_dict["polynomials"]["polynomial2"]["coeffs"]["object"]) == tuple(self.coeffs_high) diff --git a/test/rmgpy/thermo/thermodataTest.py b/test/rmgpy/thermo/thermodataTest.py index e3a633c2dac..cca336e1332 100644 --- a/test/rmgpy/thermo/thermodataTest.py +++ b/test/rmgpy/thermo/thermodataTest.py @@ -53,6 +53,7 @@ def setup_class(self): self.Tmin = 100.0 self.Tmax = 3000.0 self.E0 = -782292.0 + self.thermo_coverage_dependence = {'OX':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}} self.comment = "C2H6" self.thermodata = ThermoData( Tdata=(self.Tdata, "K"), @@ -64,6 +65,7 @@ def setup_class(self): Tmin=(self.Tmin, "K"), Tmax=(self.Tmax, "K"), E0=(self.E0, "J/mol"), + thermo_coverage_dependence=self.thermo_coverage_dependence, comment=self.comment, ) @@ -130,6 +132,12 @@ def test_comment(self): Test that the ThermoData comment property was properly set. """ assert self.thermodata.comment == self.comment + + def test_thermo_coverage_dependence(self): + """ + Test that the ThermoData thermo_coverage_dependence property was properly set. + """ + assert repr(self.thermodata.thermo_coverage_dependence) == repr(self.thermo_coverage_dependence) def test_is_temperature_valid(self): """ @@ -261,6 +269,7 @@ def test_pickle(self): assert round(abs(self.thermodata.E0.value - thermodata.E0.value), 4) == 0 assert self.thermodata.E0.units == thermodata.E0.units assert self.thermodata.label == thermodata.label + assert repr(self.thermodata.thermo_coverage_dependence) == repr(thermodata.thermo_coverage_dependence) assert self.thermodata.comment == thermodata.comment def test_repr(self): @@ -295,6 +304,7 @@ def test_repr(self): assert round(abs(self.thermodata.E0.value - thermodata.E0.value), 4) == 0 assert self.thermodata.E0.units == thermodata.E0.units assert self.thermodata.label == thermodata.label + assert repr(self.thermodata.thermo_coverage_dependence) == repr(thermodata.thermo_coverage_dependence) assert self.thermodata.comment == thermodata.comment def test_is_all_zeros(self): diff --git a/test/rmgpy/thermo/wilhoitTest.py b/test/rmgpy/thermo/wilhoitTest.py index 4b74fd72a76..15d76d34fa1 100644 --- a/test/rmgpy/thermo/wilhoitTest.py +++ b/test/rmgpy/thermo/wilhoitTest.py @@ -45,7 +45,6 @@ class TestWilhoit: """ Contains unit tests of the :class:`Wilhoit` class. """ - def setup_class(self): self.Cp0 = 4.0 self.CpInf = 21.5 @@ -59,6 +58,7 @@ def setup_class(self): self.Tmin = 300.0 self.Tmax = 3000.0 self.comment = "C2H6" + self.thermo_coverage_dependence = {'OX':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}} self.wilhoit = Wilhoit( Cp0=(self.Cp0 * constants.R, "J/(mol*K)"), CpInf=(self.CpInf * constants.R, "J/(mol*K)"), @@ -71,6 +71,7 @@ def setup_class(self): S0=(self.S0 * constants.R, "J/(mol*K)"), Tmin=(self.Tmin, "K"), Tmax=(self.Tmax, "K"), + thermo_coverage_dependence=self.thermo_coverage_dependence, comment=self.comment, ) @@ -159,6 +160,12 @@ def test_comment(self): Test that the Wilhoit comment property was properly set. """ assert self.wilhoit.comment == self.comment + + def test_thermo_coverage_dependence(self): + """ + Test that the Wilhoit thermo_coverage_dependence property was properly set. + """ + assert repr(self.wilhoit.thermo_coverage_dependence) == repr(self.thermo_coverage_dependence) def test_is_temperature_valid(self): """ @@ -287,6 +294,7 @@ def test_pickle(self): assert self.wilhoit.Tmax.units == wilhoit.Tmax.units assert round(abs(self.wilhoit.E0.value - wilhoit.E0.value), 4) == 0 assert self.wilhoit.E0.units == wilhoit.E0.units + assert repr(self.wilhoit.thermo_coverage_dependence) == repr(wilhoit.thermo_coverage_dependence) assert self.wilhoit.comment == wilhoit.comment def test_repr(self): @@ -318,6 +326,7 @@ def test_repr(self): assert self.wilhoit.Tmax.units == wilhoit.Tmax.units assert round(abs(self.wilhoit.E0.value - wilhoit.E0.value), 1) == 0 assert self.wilhoit.E0.units == wilhoit.E0.units + assert repr(self.wilhoit.thermo_coverage_dependence) == repr(wilhoit.thermo_coverage_dependence) assert self.wilhoit.comment == wilhoit.comment def test_fit_to_data(self): @@ -381,6 +390,7 @@ def test_to_wilhoit(self): spc = Species().from_smiles("CC") spc.get_thermo_data() + spc.thermo.thermo_coverage_dependence = self.thermo_coverage_dependence T = 1350.0 # not 298K! @@ -448,6 +458,16 @@ def test_wilhoit_as_dict(self): "value": 178.76114800000002, }, "class": "Wilhoit", + 'thermo_coverage_dependence': {'OX': { + 'model': 'polynomial', + 'enthalpy-coefficients': [{'class': 'ScalarQuantity', 'value': 1.0}, + {'class': 'ScalarQuantity', 'value': 2.0}, + {'class': 'ScalarQuantity', 'value': 3.0} + ], + 'entropy-coefficients': [{'class': 'ScalarQuantity', 'value': 1.0}, + {'class': 'ScalarQuantity', 'value': 2.0}, + {'class': 'ScalarQuantity', 'value': 3.0} + ]}} } def test_make_wilhoit(self): From cdf11f3a3260119df8df794ba202b39ca3370844 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Wed, 13 Mar 2024 16:11:59 -0400 Subject: [PATCH 11/49] check thermo_coverage_dependence is not empty --- test/rmgpy/thermo/convertTest.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/rmgpy/thermo/convertTest.py b/test/rmgpy/thermo/convertTest.py index d35f3de4c7f..ef2f7265e31 100644 --- a/test/rmgpy/thermo/convertTest.py +++ b/test/rmgpy/thermo/convertTest.py @@ -132,6 +132,7 @@ def test_convert_wilhoit_to_nasa(self): s_nasa = nasa.get_entropy(T) assert abs(s_nasa - s_wilhoit) < 1e0 assert abs(wilhoit.E0.value_si - nasa.E0.value_si) < 1e1 + assert repr(nasa.thermo_coverage_dependence) != {} assert repr(nasa.thermo_coverage_dependence) == repr(wilhoit.thermo_coverage_dependence) def test_convert_wilhoit_to_thermo_data(self): @@ -153,6 +154,7 @@ def test_convert_wilhoit_to_thermo_data(self): s_thermodata = thermodata.get_entropy(T) assert round(abs(s_thermodata - s_wilhoit), 4) == 0 assert abs(wilhoit.E0.value_si - thermodata.E0.value_si) < 1e1 + assert repr(thermodata.thermo_coverage_dependence) != {} assert repr(thermodata.thermo_coverage_dependence) == repr(wilhoit.thermo_coverage_dependence) def test_convert_nasa_to_wilhoit(self): @@ -173,6 +175,7 @@ def test_convert_nasa_to_wilhoit(self): s_nasa = nasa.get_entropy(T) assert abs(s_nasa - s_wilhoit) < 1e0 assert abs(nasa.E0.value_si - wilhoit.E0.value_si) < 2e1 + assert repr(wilhoit.thermo_coverage_dependence) != {} assert repr(wilhoit.thermo_coverage_dependence) == repr(nasa.thermo_coverage_dependence) def test_convert_nasa_to_thermo_data(self): @@ -194,6 +197,7 @@ def test_convert_nasa_to_thermo_data(self): s_nasa = nasa.get_entropy(T) assert round(abs(s_nasa - s_thermodata), 4) == 0 assert abs(nasa.E0.value_si - thermodata.E0.value_si) < 1e1 + assert repr(thermodata.thermo_coverage_dependence) != {} assert repr(thermodata.thermo_coverage_dependence) == repr(nasa.thermo_coverage_dependence) def test_convert_thermo_data_to_wilhoit(self): @@ -215,6 +219,7 @@ def test_convert_thermo_data_to_wilhoit(self): s_thermodata = thermodata.get_entropy(T) assert round(abs(s_thermodata - s_wilhoit), 3) == 0 assert abs(thermodata.E0.value_si - wilhoit.E0.value_si) < 1e1 + assert repr(wilhoit.thermo_coverage_dependence) != {} assert repr(wilhoit.thermo_coverage_dependence) == repr(thermodata.thermo_coverage_dependence) def test_convert_thermo_data_to_nasa(self): @@ -236,6 +241,7 @@ def test_convert_thermo_data_to_nasa(self): s_nasa = nasa.get_entropy(T) assert abs(s_nasa - s_thermodata) < 1e0 assert abs(thermodata.E0.value_si - nasa.E0.value_si) < 1e1 + assert repr(nasa.thermo_coverage_dependence) != {} assert repr(nasa.thermo_coverage_dependence) == repr(thermodata.thermo_coverage_dependence) def test_wilhoit_nasa_wilhoit(self): @@ -257,6 +263,7 @@ def test_wilhoit_nasa_wilhoit(self): s_2 = wilhoit2.get_entropy(T) assert abs(s_1 - s_2) < 1e0 assert abs(wilhoit1.E0.value_si - wilhoit2.E0.value_si) < 1e1 + assert repr(wilhoit1.thermo_coverage_dependence) != {} assert repr(wilhoit1.thermo_coverage_dependence) == repr(wilhoit2.thermo_coverage_dependence) def test_wilhoit_thermo_data_wilhoit(self): @@ -278,4 +285,5 @@ def test_wilhoit_thermo_data_wilhoit(self): s_2 = wilhoit2.get_entropy(T) assert abs(s_1 - s_2) < 1e0 assert abs(wilhoit1.E0.value_si - wilhoit2.E0.value_si) < 1e1 + assert repr(wilhoit1.thermo_coverage_dependence) != {} assert repr(wilhoit1.thermo_coverage_dependence) == repr(wilhoit2.thermo_coverage_dependence) From d927e5adac56a98919d0c8031edfc84cc8b7a104 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Wed, 13 Mar 2024 19:19:42 -0400 Subject: [PATCH 12/49] Fix the bug for calculating thermo coverage dependent correct values --- rmgpy/solver/surface.pyx | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index 54151c49938..4b15a0240f8 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -205,9 +205,7 @@ cdef class SurfaceReactor(ReactionSystem): 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: - # think about the data structure of thermo coverage dependence - # should not be a list, a dictionary would be better - for spec, parameters in sp.thermo.coverage_dependence.items(): + for spec, parameters in sp.thermo.thermo_coverage_dependence.items(): species_index = self.species_index[spec] try: list_of_thermo_coverage_deps = self.thermo_coverage_dependencies[species_index] @@ -216,6 +214,10 @@ cdef class SurfaceReactor(ReactionSystem): self.thermo_coverage_dependencies[sp_index] = list_of_thermo_coverage_deps # need to specify the entropy and enthalpy models # linear, piecewise linear, polynomial, interpolative models + if parameters['model'] == "polynomial": + # for the case of polynomial, we need to insert a 0 for the constant term + parameters['enthalpy-coefficients'] = [0]+[x.value_si for x in parameters['enthalpy-coefficients']] + parameters['entropy-coefficients'] = [0]+[x.value_si for x in parameters['entropy-coefficients']] list_of_thermo_coverage_deps.append((sp_index, parameters)) """ @@ -478,14 +480,14 @@ cdef class SurfaceReactor(ReactionSystem): if parameters['model'] == "linear": pass elif parameters['model'] == "polynomial": - enthalpy_cov_correction = np.polynomial.polynomial.polyval(surface_site_fraction, parameters['enthalpy-coefficients'].insert(0,0)) # insert 0 for the constant term - entropy_cov_correction = np.polynomial.polynomial.polyval(surface_site_fraction, parameters['entropy-coefficients'].insert(0,0)) + enthalpy_cov_correction = np.polynomial.polynomial.polyval(surface_site_fraction, parameters['enthalpy-coefficients']) # insert 0 for the constant term + entropy_cov_correction = np.polynomial.polynomial.polyval(surface_site_fraction, parameters['entropy-coefficients']) free_energy_coverage_corrections[j] += enthalpy_cov_correction - self.T.value_si * entropy_cov_correction elif parameters['model'] == "piecewise-linear": pass elif parameters['model'] == "interpolative": pass - corrected_K_eq = copy.deepcopy(self.K_eq) + corrected_K_eq = copy.deepcopy(self.Keq) # correct the K_eq for j in range(ir.shape[0]): if ir[j, 0] >= num_core_species or ir[j, 1] >= num_core_species or ir[j, 2] >= num_core_species: From 060f252e73b69d7cb02c1bedbcc78c687d1791e1 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Fri, 15 Mar 2024 17:58:21 -0400 Subject: [PATCH 13/49] add a big matrix for species Gibbs free energy corrections --- rmgpy/solver/surface.pyx | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index 4b15a0240f8..ab334546b82 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -67,6 +67,7 @@ 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 bint coverage_dependence cdef public dict coverage_dependencies @@ -173,6 +174,8 @@ 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), 4), dtype=np.float64) + 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) @@ -206,7 +209,10 @@ cdef class SurfaceReactor(ReactionSystem): 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(): - species_index = self.species_index[spec] + try: + species_index = self.species_index[spec] + except KeyError: + logging.warning("Species {} is not in the species list yet, skip the thermodynamic coverage effect estimation!".format(spec)) try: list_of_thermo_coverage_deps = self.thermo_coverage_dependencies[species_index] except KeyError: # doesn't exist yet @@ -413,8 +419,6 @@ cdef class SurfaceReactor(ReactionSystem): cdef list list_of_coverage_deps cdef double surface_site_fraction, total_sites, a, m, E - - ir = self.reactant_indices ip = self.product_indices equilibrium_constants = self.Keq @@ -500,12 +504,12 @@ cdef class SurfaceReactor(ReactionSystem): corrected_K_eq[j] *= np.exp((free_energy_coverage_corrections[ir[j, 0]] + free_energy_coverage_corrections[ir[j, 1]] + free_energy_coverage_corrections[ir[j, 2]]) / (constants.R * self.T.value_si)) if ip[j, 0] >= num_core_species or ip[j, 1] >= num_core_species or ip[j, 2] >= num_core_species: pass - elif ip[j, 1] == -1: # only one reactant - corrected_K_eq[j] /= np.exp(free_energy_coverage_corrections[ir[j, 0]] / (constants.R * self.T.value_si)) - elif ip[j, 2] == -1: # only two reactants - corrected_K_eq[j] /= np.exp((free_energy_coverage_corrections[ir[j, 0]] + free_energy_coverage_corrections[ir[j, 1]]) / (constants.R * self.T.value_si)) - else: # three reactants!! (really?) - corrected_K_eq[j] /= np.exp((free_energy_coverage_corrections[ir[j, 0]] + free_energy_coverage_corrections[ir[j, 1]] + free_energy_coverage_corrections[ir[j, 2]]) / (constants.R * self.T.value_si)) + elif ip[j, 1] == -1: # only one product + corrected_K_eq[j] /= np.exp(free_energy_coverage_corrections[ip[j, 0]] / (constants.R * self.T.value_si)) + elif ip[j, 2] == -1: # only two products + corrected_K_eq[j] /= np.exp((free_energy_coverage_corrections[ip[j, 0]] + free_energy_coverage_corrections[ip[j, 1]]) / (constants.R * self.T.value_si)) + else: # three products!! (really?) + corrected_K_eq[j] /= np.exp((free_energy_coverage_corrections[ip[j, 0]] + free_energy_coverage_corrections[ip[j, 1]] + free_energy_coverage_corrections[ip[j, 2]]) / (constants.R * self.T.value_si)) kr = kf / corrected_K_eq # Coverage dependence coverage_corrections = np.ones_like(kf, float) From 518859e31b33e91118b61ecafe483af22bbfb8ef Mon Sep 17 00:00:00 2001 From: 12Chao Date: Fri, 15 Mar 2024 21:25:26 -0400 Subject: [PATCH 14/49] correct Gibbs free energy for each species based on its 3-parameter polynomial thermo coverage depedent model --- rmgpy/solver/surface.pyx | 72 +++++++++++----------------------------- 1 file changed, 19 insertions(+), 53 deletions(-) diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index ab334546b82..b943df82636 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -174,7 +174,7 @@ 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), 4), dtype=np.float64) + cdef np.ndarray thermo_coeff_matrix = np.zeros((len(self.species_index), len(self.species_index), 6), dtype=np.float64) 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) @@ -211,30 +211,11 @@ cdef class SurfaceReactor(ReactionSystem): for spec, parameters in sp.thermo.thermo_coverage_dependence.items(): try: species_index = self.species_index[spec] + thermo_polynomials = parameters['enthalpy-coefficients'] + parameters['entropy-coefficients'] + self.thermo_coeff_matrix[sp_index, species_index] = [x.value_si for x in thermo_polynomials] except KeyError: logging.warning("Species {} is not in the species list yet, skip the thermodynamic coverage effect estimation!".format(spec)) - try: - list_of_thermo_coverage_deps = self.thermo_coverage_dependencies[species_index] - except KeyError: # doesn't exist yet - list_of_thermo_coverage_deps = [] - self.thermo_coverage_dependencies[sp_index] = list_of_thermo_coverage_deps - # need to specify the entropy and enthalpy models - # linear, piecewise linear, polynomial, interpolative models - if parameters['model'] == "polynomial": - # for the case of polynomial, we need to insert a 0 for the constant term - parameters['enthalpy-coefficients'] = [0]+[x.value_si for x in parameters['enthalpy-coefficients']] - parameters['entropy-coefficients'] = [0]+[x.value_si for x in parameters['entropy-coefficients']] - list_of_thermo_coverage_deps.append((sp_index, parameters)) - - """ - self.thermo_coverage_dependencies[2] = [(3, {"model":"linear","enthalpy-coefficients":[], "entropy-coefficients":[]}),] - self.thermo_coverage_dependencies[2] = [(3, {"model":"polynomial","enthalpy-coefficients":[], "entropy-coefficients":[]}),] - self.thermo_coverage_dependencies[2] = [(3, {"model":"piecewise-linear","enthalpy-coefficients":{"enthalpy_high":, "enthalpy_low":, "enthalpy_change":}, - "entropy-coefficients":{"entropy_high":, "entropy_low":, "entropy_change"}, "Cp":{"heat-capacity-a":, "heat-capacity-b":}}),] - self.thermo_coverage_dependencies[2] = {(3, {"model":"interpolative","enthalpy-coefficients":{"enthalpy-coverages":, "enthalpies":}, "entropy-coefficients":{"entropy-coverages":, "entropies":}}),} - means that Species with index 2 in the current simulation is used in - Species 3 with parameters for linear, polynomial, piecewise-linear, and interpolative models - """ + self.species_on_surface = species_on_surface self.reactions_on_surface = reactions_on_surface @@ -418,7 +399,6 @@ 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 - ir = self.reactant_indices ip = self.product_indices equilibrium_constants = self.Keq @@ -452,7 +432,6 @@ 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 @@ -462,35 +441,22 @@ cdef class SurfaceReactor(ReactionSystem): core_species_concentrations[j] = C[j] # Thermodynamic coverage dependence - free_energy_coverage_corrections = np.zeros(len(self.species_index), float) # length of core + edge species if self.thermo_coverage_dependence: - """ - self.thermo_coverage_dependencies[2] = [(3, {"model":"linear","enthalpy-coefficients":[], "entropy-coefficients":[]}),] - self.thermo_coverage_dependencies[2] = [(3, {"model":"polynomial","enthalpy-coefficients":[], "entropy-coefficients":[]}),] - self.thermo_coverage_dependencies[2] = [(3, {"model":"piecewise-linear","enthalpy-coefficients":{"enthalpy_high":, "enthalpy_low":, "enthalpy_change":}, - "entropy-coefficients":{"entropy_high":, "entropy_low":, "entropy_change"}, "Cp":{"heat-capacity-a":, "heat-capacity-b":}}),] - self.thermo_coverage_dependencies[2] = {(3, {"model":"interpolative","enthalpy-coefficients":{"enthalpy-coverages":, "enthalpies":}, "entropy-coefficients":{"entropy-coverages":, "entropies":}}),} - means that Species with index 2 in the current simulation is used in - Species 3 with parameters for linear, polynomial, piecewise-linear, and interpolative models - """ - for i, list_of_thermo_coverage_deps in self.thermo_coverage_dependencies.items(): - surface_site_fraction = N[i] / total_sites - if surface_site_fraction < 1e-15: - continue - for j, parameters in list_of_thermo_coverage_deps: - # Species i, Species j - # need to specify the entropy and enthalpy models - # linear, piecewise linear, polynomial, interpolative models - if parameters['model'] == "linear": - pass - elif parameters['model'] == "polynomial": - enthalpy_cov_correction = np.polynomial.polynomial.polyval(surface_site_fraction, parameters['enthalpy-coefficients']) # insert 0 for the constant term - entropy_cov_correction = np.polynomial.polynomial.polyval(surface_site_fraction, parameters['entropy-coefficients']) - free_energy_coverage_corrections[j] += enthalpy_cov_correction - self.T.value_si * entropy_cov_correction - elif parameters['model'] == "piecewise-linear": - pass - elif parameters['model'] == "interpolative": - pass + coverages = [] + for i in range(len(N)): + if species_on_surface[i]: + surface_site_fraction = N[i] / total_sites + else: + surface_site_fraction = 0 + coverages.append(surface_site_fraction) + coverages = np.array(coverages) + thermo_dep_coverage = np.stack([coverages, coverages**2, coverages**3, -self.T.value_si*coverages, -self.T.value_si*coverages**2, -self.T.value_si*coverages**3]) + free_energy_coverage_corrections = [] + for matrix in self.thermo_coeff_matrix: + sp_free_energy_correction = np.diag(np.dot(matrix, thermo_dep_coverage)).sum() + free_energy_coverage_corrections.append(sp_free_energy_correction) + free_energy_coverage_corrections = np.array(free_energy_coverage_corrections) + corrected_K_eq = copy.deepcopy(self.Keq) # correct the K_eq for j in range(ir.shape[0]): From e29275553d439a0e04e0957b6c35820d25441680 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 18 Mar 2024 19:24:25 -0400 Subject: [PATCH 15/49] use stoichiometric matrix to calculate the correct value for Keq --- rmgpy/solver/surface.pyx | 57 +++++++++++++++++++++++++--------------- 1 file changed, 36 insertions(+), 21 deletions(-) diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index b943df82636..cddd981d274 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -68,6 +68,7 @@ cdef class SurfaceReactor(ReactionSystem): 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 @@ -175,7 +176,9 @@ 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) self.thermo_coeff_matrix = thermo_coeff_matrix + self.stoi_matrix = stoi_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) @@ -215,7 +218,36 @@ cdef class SurfaceReactor(ReactionSystem): self.thermo_coeff_matrix[sp_index, species_index] = [x.value_si for x in thermo_polynomials] except KeyError: logging.warning("Species {} is not in the species list yet, skip the thermodynamic coverage effect estimation!".format(spec)) - + + 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 @@ -455,28 +487,11 @@ cdef class SurfaceReactor(ReactionSystem): for matrix in self.thermo_coeff_matrix: sp_free_energy_correction = np.diag(np.dot(matrix, thermo_dep_coverage)).sum() free_energy_coverage_corrections.append(sp_free_energy_correction) - free_energy_coverage_corrections = np.array(free_energy_coverage_corrections) - + rxns_free_energy_change = np.diag(np.dot(self.stoi_matrix, np.transpose(np.array([free_energy_coverage_corrections])))) corrected_K_eq = copy.deepcopy(self.Keq) - # correct the K_eq - for j in range(ir.shape[0]): - if ir[j, 0] >= num_core_species or ir[j, 1] >= num_core_species or ir[j, 2] >= num_core_species: - pass - elif ir[j, 1] == -1: # only one reactant - corrected_K_eq[j] *= np.exp(free_energy_coverage_corrections[ir[j, 0]] / (constants.R * self.T.value_si)) - elif ir[j, 2] == -1: # only two reactants - corrected_K_eq[j] *= np.exp((free_energy_coverage_corrections[ir[j, 0]] + free_energy_coverage_corrections[ir[j, 1]]) / (constants.R * self.T.value_si)) - else: # three reactants!! (really?) - corrected_K_eq[j] *= np.exp((free_energy_coverage_corrections[ir[j, 0]] + free_energy_coverage_corrections[ir[j, 1]] + free_energy_coverage_corrections[ir[j, 2]]) / (constants.R * self.T.value_si)) - if ip[j, 0] >= num_core_species or ip[j, 1] >= num_core_species or ip[j, 2] >= num_core_species: - pass - elif ip[j, 1] == -1: # only one product - corrected_K_eq[j] /= np.exp(free_energy_coverage_corrections[ip[j, 0]] / (constants.R * self.T.value_si)) - elif ip[j, 2] == -1: # only two products - corrected_K_eq[j] /= np.exp((free_energy_coverage_corrections[ip[j, 0]] + free_energy_coverage_corrections[ip[j, 1]]) / (constants.R * self.T.value_si)) - else: # three products!! (really?) - corrected_K_eq[j] /= np.exp((free_energy_coverage_corrections[ip[j, 0]] + free_energy_coverage_corrections[ip[j, 1]] + free_energy_coverage_corrections[ip[j, 2]]) / (constants.R * self.T.value_si)) + 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) if self.coverage_dependence: From 64b6e20a9720be62c20065f6b4fd5bf0c55e7dd8 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Wed, 20 Mar 2024 12:05:02 -0400 Subject: [PATCH 16/49] Add unit tests for thermo coverage dependent models in surface solver tests --- test/rmgpy/solver/solverSurfaceTest.py | 394 +++++++++++++++++++++++++ 1 file changed, 394 insertions(+) diff --git a/test/rmgpy/solver/solverSurfaceTest.py b/test/rmgpy/solver/solverSurfaceTest.py index 36c08b4ad83..dc9ee25b3fa 100644 --- a/test/rmgpy/solver/solverSurfaceTest.py +++ b/test/rmgpy/solver/solverSurfaceTest.py @@ -762,3 +762,397 @@ def test_solve_ch3_coverage_dependence(self): # Check that coverages are different assert not np.allclose(y, y_off) assert not np.allclose(species_rates, species_rates_off) + + def test_solve_h2_thermo_coverage_dependence(self): + """ + Test the surface batch reactor can properly apply thermo coverage dependent parameters + with the dissociative adsorption of H2. + + Here we choose a kinetic model consisting of the dissociative adsorption reaction + H2 + 2X <=> 2 HX + We use a SurfaceArrhenius for the rate expression. + """ + h2 = Species( + molecule=[Molecule().from_smiles("[H][H]")], + thermo=ThermoData( + Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"), + Cpdata=( + [6.955, 6.955, 6.956, 6.961, 7.003, 7.103, 7.502], + "cal/(mol*K)", + ), + H298=(0, "kcal/mol"), + S298=(31.129, "cal/(mol*K)"), + ), + ) + + x = Species( + molecule=[Molecule().from_adjacency_list("1 X u0 p0")], + thermo=ThermoData( + Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"), + Cpdata=([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], "cal/(mol*K)"), + H298=(0.0, "kcal/mol"), + S298=(0.0, "cal/(mol*K)"), + ), + ) + + hx = Species( + molecule=[Molecule().from_adjacency_list("1 H u0 p0 {2,S} \n 2 X u0 p0 {1,S}")], + thermo=ThermoData( + Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"), + Cpdata=([1.50, 2.58, 3.40, 4.00, 4.73, 5.13, 5.57], "cal/(mol*K)"), + H298=(-11.26, "kcal/mol"), + S298=(0.44, "cal/(mol*K)"), + ), + ) + hx.thermo.thermo_coverage_dependence = {hx:{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,5,3]},} + + rxn1 = Reaction( + reactants=[h2, x, x], + products=[hx, hx], + kinetics=SurfaceArrhenius( + A=(9.05e18, "cm^5/(mol^2*s)"), + n=0.5, + Ea=(5.0, "kJ/mol"), + T0=(1.0, "K"), + ), + ) + + rxn2 = Reaction( + reactants=[h2, x, x], + products=[hx, hx], + kinetics=SurfaceArrhenius( + A=(9.05e-18, "cm^5/(mol^2*s)"), # 1e36 times slower + n=0.5, + Ea=(5.0, "kJ/mol"), + T0=(1.0, "K"), + ), + ) + + core_species = [h2, x, hx] + edge_species = [] + core_reactions = [rxn1] + edge_reactions = [] + + T = 600 + P_initial = 1.0e5 + rxn_system = SurfaceReactor( + T, + P_initial, + n_sims=1, + initial_gas_mole_fractions={h2: 1.0}, + initial_surface_coverages={x: 1.0}, + surface_volume_ratio=(1e1, "m^-1"), + surface_site_density=(2.72e-9, "mol/cm^2"), + coverage_dependence=True, + termination=[], + ) + + rxn_system.initialize_model(core_species, core_reactions, edge_species, edge_reactions) + + tlist = np.logspace(-13, -5, 81, dtype=float) + + assert isinstance(hx.thermo.thermo_coverage_dependence, dict) # check to make sure coverage_dependence is still the correct type + for species, parameters in hx.thermo.thermo_coverage_dependence.items(): + assert isinstance(species, Species) # species should be a Species + assert isinstance(parameters, dict) + assert parameters["model"] is not None + assert parameters["enthalpy-coefficients"] is not None + assert parameters["entropy-coefficients"] is not None + # Integrate to get the solution at each time point + t = [] + y = [] + reaction_rates = [] + species_rates = [] + start_time = time.time() + for t1 in tlist: + rxn_system.advance(t1) + t.append(rxn_system.t) + # You must make a copy of y because it is overwritten by DASSL at + # each call to advance() + y.append(rxn_system.y.copy()) + reaction_rates.append(rxn_system.core_reaction_rates.copy()) + species_rates.append(rxn_system.core_species_rates.copy()) + run_time = time.time() - start_time + print(f"Simulation took {run_time:.3e} seconds") + + # Convert the solution vectors to np arrays + t = np.array(t, float) + y = np.array(y, float) + reaction_rates = np.array(reaction_rates, float) + species_rates = np.array(species_rates, float) + total_sites = y[0, 1] + + # Check that we're computing the species fluxes correctly + for i in range(t.shape[0]): + assert abs(reaction_rates[i, 0] - -1.0 * species_rates[i, 0]) < abs( + 1e-6 * reaction_rates[i, 0] + ) + assert abs(reaction_rates[i, 0] - -0.5 * species_rates[i, 1]) < abs( + 1e-6 * reaction_rates[i, 0] + ) + assert abs(reaction_rates[i, 0] - 0.5 * species_rates[i, 2]) < abs( + 1e-6 * reaction_rates[i, 0] + ) + + # Check that we've reached equilibrium + assert abs(reaction_rates[-1, 0] - 0.0) < 1e-2 + + def test_solve_ch3_thermo_coverage_dependence(self): + """ + Test the surface batch reactor can properly apply coverage dependent parameters + with the nondissociative adsorption of CH3 + + Here we choose a kinetic model consisting of the adsorption reaction + CH3 + X <=> CH3X + We use a sticking coefficient for the rate expression. + """ + + ch3 = Species( + molecule=[Molecule().from_smiles("[CH3]")], + thermo=NASA( + polynomials=[ + NASAPolynomial( + coeffs=[ + 3.91547, + 0.00184155, + 3.48741e-06, + -3.32746e-09, + 8.49953e-13, + 16285.6, + 0.351743, + ], + Tmin=(100, "K"), + Tmax=(1337.63, "K"), + ), + NASAPolynomial( + coeffs=[ + 3.54146, + 0.00476786, + -1.82148e-06, + 3.28876e-10, + -2.22545e-14, + 16224, + 1.66032, + ], + Tmin=(1337.63, "K"), + Tmax=(5000, "K"), + ), + ], + Tmin=(100, "K"), + Tmax=(5000, "K"), + E0=(135.382, "kJ/mol"), + comment="""Thermo library: primaryThermoLibrary + radical(CH3)""", + ), + molecular_weight=(15.0345, "amu"), + ) + + x = Species( + molecule=[Molecule().from_adjacency_list("1 X u0 p0")], + thermo=NASA( + polynomials=[ + NASAPolynomial(coeffs=[0, 0, 0, 0, 0, 0, 0], Tmin=(298, "K"), Tmax=(1000, "K")), + NASAPolynomial(coeffs=[0, 0, 0, 0, 0, 0, 0], Tmin=(1000, "K"), Tmax=(2000, "K")), + ], + Tmin=(298, "K"), + Tmax=(2000, "K"), + E0=(-6.19426, "kJ/mol"), + comment="""Thermo library: surfaceThermo""", + ), + ) + + ch3x = Species( + molecule=[ + Molecule().from_adjacency_list( + """1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S} +2 H u0 p0 c0 {1,S} +3 H u0 p0 c0 {1,S} +4 H u0 p0 c0 {1,S} +5 X u0 p0 c0 {1,S}""" + ) + ], + thermo=NASA( + polynomials=[ + NASAPolynomial( + coeffs=[ + -0.552219, + 0.026442, + -3.55617e-05, + 2.60044e-08, + -7.52707e-12, + -4433.47, + 0.692144, + ], + Tmin=(298, "K"), + Tmax=(1000, "K"), + ), + NASAPolynomial( + coeffs=[ + 3.62557, + 0.00739512, + -2.43797e-06, + 1.86159e-10, + 3.6485e-14, + -5187.22, + -18.9668, + ], + Tmin=(1000, "K"), + Tmax=(2000, "K"), + ), + ], + Tmin=(298, "K"), + Tmax=(2000, "K"), + E0=(-39.1285, "kJ/mol"), + comment="""Thermo library: surfaceThermoNi111""", + ), + ) + ch3x.thermo.thermo_coverage_dependence = {ch3x:{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,5,3]},} + + rxn1 = Reaction( + reactants=[ch3, x], + products=[ch3x], + kinetics=StickingCoefficient( + A=0.1, + n=0, + Ea=(0, "kcal/mol"), + T0=(1, "K"), + Tmin=(200, "K"), + Tmax=(3000, "K"), + coverage_dependence={x: {"a": 0.0, "m": -1.0, "E": (0.0, "J/mol")}}, + comment="""Exact match found for rate rule (Adsorbate;VacantSite)""", + ), + ) + core_species = [ch3, x, ch3x] + edge_species = [] + core_reactions = [rxn1] + edge_reactions = [] + + T = 800.0 + P_initial = 1.0e5 + rxn_system = SurfaceReactor( + T, + P_initial, + n_sims=1, + initial_gas_mole_fractions={ch3: 1.0}, + initial_surface_coverages={x: 1.0}, + surface_volume_ratio=(1.0, "m^-1"), + surface_site_density=(2.72e-9, "mol/cm^2"), + coverage_dependence=True, + termination=[], + ) + # in chemkin, the sites are mostly occupied in about 1e-8 seconds. + + rxn_system.initialize_model(core_species, core_reactions, edge_species, edge_reactions) + + tlist = np.logspace(-13, -5, 81, dtype=float) + + print("Surface site density:", rxn_system.surface_site_density.value_si) + + print( + "rxn1 rate coefficient", + rxn1.get_surface_rate_coefficient(rxn_system.T.value_si, rxn_system.surface_site_density.value_si), + ) + + assert isinstance(ch3.thermo.thermo_coverage_dependence, dict) # check to make sure coverage_dependence is still the correct type + for species, parameters in ch3.thermo.thermo_coverage_dependence.items(): + assert isinstance(species, Species) # species should be a Species + assert isinstance(parameters, dict) + assert parameters["model"] is not None + assert parameters["enthalpy-coefficients"] is not None + assert parameters["entropy-coefficients"] is not None + + # Integrate to get the solution at each time point + t = [] + y = [] + reaction_rates = [] + species_rates = [] + t.append(rxn_system.t) + # You must make a copy of y because it is overwritten by DASSL at + # each call to advance() + y.append(rxn_system.y.copy()) + reaction_rates.append(rxn_system.core_reaction_rates.copy()) + species_rates.append(rxn_system.core_species_rates.copy()) + print("time: ", t) + print("moles:", y) + print("reaction rates:", reaction_rates) + print("species rates:", species_rates) + start_time = time.time() + for t1 in tlist: + rxn_system.advance(t1) + t.append(rxn_system.t) + # You must make a copy of y because it is overwritten by DASSL at + # each call to advance() + y.append(rxn_system.y.copy()) + reaction_rates.append(rxn_system.core_reaction_rates.copy()) + species_rates.append(rxn_system.core_species_rates.copy()) + run_time = time.time() - start_time + print(f"Simulation took {run_time:.3e} seconds") + + # Convert the solution vectors to np arrays + t = np.array(t, float) + y = np.array(y, float) + reaction_rates = np.array(reaction_rates, float) + species_rates = np.array(species_rates, float) + V = constants.R * rxn_system.T.value_si * np.sum(y) / rxn_system.P_initial.value_si + + # Check that we're computing the species fluxes correctly + for i in range(t.shape[0]): + assert abs(reaction_rates[i, 0] - -species_rates[i, 0]) < abs( + 1e-6 * reaction_rates[i, 0] + ) + assert abs(reaction_rates[i, 0] - -species_rates[i, 1]) < abs( + 1e-6 * reaction_rates[i, 0] + ) + assert abs(reaction_rates[i, 0] - species_rates[i, 2]) < abs( + 1e-6 * reaction_rates[i, 0] + ) + + # Check that we've reached equilibrium by the end + assert abs(reaction_rates[-1, 0] - 0.0) < 1e-2 + + # Run model with Covdep off so we can test that it is actually being implemented + rxn_system = SurfaceReactor( + T, + P_initial, + n_sims=1, + initial_gas_mole_fractions={ch3: 1.0}, + initial_surface_coverages={x: 1.0}, + surface_volume_ratio=(1.0, "m^-1"), + surface_site_density=(2.72e-9, "mol/cm^2"), + termination=[], + ) + + rxn_system.initialize_model(core_species, core_reactions, edge_species, edge_reactions) + + tlist = np.logspace(-13, -5, 81, dtype=float) + + # Integrate to get the solution at each time point + t = [] + y_off = [] + species_rates_off = [] + t.append(rxn_system.t) + + # You must make a copy of y because it is overwritten by DASSL at + # each call to advance() + y_off.append(rxn_system.y.copy()) + species_rates_off.append(rxn_system.core_species_rates.copy()) + for t1 in tlist: + rxn_system.advance(t1) + t.append(rxn_system.t) + # You must make a copy of y because it is overwritten by DASSL at + # each call to advance() + y_off.append(rxn_system.y.copy()) + species_rates_off.append(rxn_system.core_species_rates.copy()) + run_time = time.time() - start_time + print(f"Simulation took {run_time:.3e} seconds") + + # Convert the solution vectors to np arrays + t = np.array(t, float) + y_off = np.array(y_off, float) + species_rates_off = np.array(species_rates_off, float) + + # Check that we've reached equilibrium + assert abs(species_rates_off[-1, 0] - 0.0) < 1e-2 + + # Check that coverages are different + assert not np.allclose(y, y_off) + assert not np.allclose(species_rates, species_rates_off) From e81f9472d98be11aefd745bd6af1c1e3abc22f3f Mon Sep 17 00:00:00 2001 From: 12Chao Date: Sun, 24 Mar 2024 19:59:55 -0400 Subject: [PATCH 17/49] use isomorphic check to check if thermo coverage dependent species are in the model --- rmgpy/solver/surface.pyx | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index cddd981d274..20b01311bf3 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -44,6 +44,7 @@ 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): """ @@ -178,7 +179,6 @@ cdef class SurfaceReactor(ReactionSystem): 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) self.thermo_coeff_matrix = thermo_coeff_matrix - self.stoi_matrix = stoi_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) @@ -212,13 +212,14 @@ cdef class SurfaceReactor(ReactionSystem): 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(): - try: - species_index = self.species_index[spec] - thermo_polynomials = parameters['enthalpy-coefficients'] + parameters['entropy-coefficients'] - self.thermo_coeff_matrix[sp_index, species_index] = [x.value_si for x in thermo_polynomials] - except KeyError: - logging.warning("Species {} is not in the species list yet, skip the thermodynamic coverage effect estimation!".format(spec)) - + 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] + thermo_polynomials = parameters['enthalpy-coefficients'] + parameters['entropy-coefficients'] + self.thermo_coeff_matrix[sp_index, species_index] = [x.value_si 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 From 5c9f83aaf52076144a2326a7b24d192b7d8fb8ad Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 25 Mar 2024 12:59:18 -0400 Subject: [PATCH 18/49] delete unused instances --- rmgpy/solver/surface.pyx | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index 20b01311bf3..becdd844f46 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -74,7 +74,6 @@ cdef class SurfaceReactor(ReactionSystem): cdef public bint coverage_dependence cdef public dict coverage_dependencies cdef public bint thermo_coverage_dependence - cdef public dict thermo_coverage_dependencies @@ -118,7 +117,6 @@ cdef class SurfaceReactor(ReactionSystem): self.n_sims = n_sims self.coverage_dependencies = {} - self.thermo_coverage_dependencies = {} def convert_initial_keys_to_species_objects(self, species_dict): """ @@ -178,7 +176,8 @@ cdef class SurfaceReactor(ReactionSystem): 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) - self.thermo_coeff_matrix = thermo_coeff_matrix + 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) From cd59089911f87a232c0f08341fd7b0e29dd68a3a Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 25 Mar 2024 13:00:30 -0400 Subject: [PATCH 19/49] Add solver unit test for thermodynamic coverage dependent models --- test/rmgpy/solver/solverSurfaceTest.py | 39 +++++++++++++++++++------- 1 file changed, 29 insertions(+), 10 deletions(-) diff --git a/test/rmgpy/solver/solverSurfaceTest.py b/test/rmgpy/solver/solverSurfaceTest.py index dc9ee25b3fa..42fb9c81404 100644 --- a/test/rmgpy/solver/solverSurfaceTest.py +++ b/test/rmgpy/solver/solverSurfaceTest.py @@ -765,7 +765,7 @@ def test_solve_ch3_coverage_dependence(self): def test_solve_h2_thermo_coverage_dependence(self): """ - Test the surface batch reactor can properly apply thermo coverage dependent parameters + Test the surface batch reactor can properly apply thermodynamic coverage dependent parameters with the dissociative adsorption of H2. Here we choose a kinetic model consisting of the dissociative adsorption reaction @@ -802,9 +802,9 @@ def test_solve_h2_thermo_coverage_dependence(self): Cpdata=([1.50, 2.58, 3.40, 4.00, 4.73, 5.13, 5.57], "cal/(mol*K)"), H298=(-11.26, "kcal/mol"), S298=(0.44, "cal/(mol*K)"), + thermo_coverage_dependence={"1 H u0 p0 {2,S} \n 2 X u0 p0 {1,S}":{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,5,3]},} ), ) - hx.thermo.thermo_coverage_dependence = {hx:{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,5,3]},} rxn1 = Reaction( reactants=[h2, x, x], @@ -844,6 +844,7 @@ def test_solve_h2_thermo_coverage_dependence(self): surface_volume_ratio=(1e1, "m^-1"), surface_site_density=(2.72e-9, "mol/cm^2"), coverage_dependence=True, + thermo_coverage_dependence=True, termination=[], ) @@ -853,11 +854,15 @@ def test_solve_h2_thermo_coverage_dependence(self): assert isinstance(hx.thermo.thermo_coverage_dependence, dict) # check to make sure coverage_dependence is still the correct type for species, parameters in hx.thermo.thermo_coverage_dependence.items(): - assert isinstance(species, Species) # species should be a Species + assert isinstance(species, str) # species should be an ajacency list assert isinstance(parameters, dict) assert parameters["model"] is not None assert parameters["enthalpy-coefficients"] is not None assert parameters["entropy-coefficients"] is not None + assert np.array_equal(rxn_system.stoi_matrix, np.array([[-1., -2., 2.]])) + thermo_coeffs = np.array([np.zeros((3,6))]*3) + thermo_coeffs[-1][-1] = [1., 2., 3., 1., 5., 3.] + assert np.array_equal(rxn_system.thermo_coeff_matrix, thermo_coeffs) # Integrate to get the solution at each time point t = [] y = [] @@ -899,7 +904,7 @@ def test_solve_h2_thermo_coverage_dependence(self): def test_solve_ch3_thermo_coverage_dependence(self): """ - Test the surface batch reactor can properly apply coverage dependent parameters + Test the surface batch reactor can properly apply thermodynamic coverage dependent parameters with the nondissociative adsorption of CH3 Here we choose a kinetic model consisting of the adsorption reaction @@ -1002,10 +1007,11 @@ def test_solve_ch3_thermo_coverage_dependence(self): Tmin=(298, "K"), Tmax=(2000, "K"), E0=(-39.1285, "kJ/mol"), + thermo_coverage_dependence={"1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S} \n 2 H u0 p0 c0 {1,S} \n 3 H u0 p0 c0 {1,S} \n 4 H u0 p0 c0 {1,S} \n 5 X u0 p0 c0 {1,S}": + {'model':'polynomial', 'enthalpy-coefficients':[1e5,2,3], "entropy-coefficients":[1,5,3]},}, comment="""Thermo library: surfaceThermoNi111""", ), ) - ch3x.thermo.thermo_coverage_dependence = {ch3x:{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,5,3]},} rxn1 = Reaction( reactants=[ch3, x], @@ -1017,7 +1023,6 @@ def test_solve_ch3_thermo_coverage_dependence(self): T0=(1, "K"), Tmin=(200, "K"), Tmax=(3000, "K"), - coverage_dependence={x: {"a": 0.0, "m": -1.0, "E": (0.0, "J/mol")}}, comment="""Exact match found for rate rule (Adsorbate;VacantSite)""", ), ) @@ -1036,7 +1041,7 @@ def test_solve_ch3_thermo_coverage_dependence(self): initial_surface_coverages={x: 1.0}, surface_volume_ratio=(1.0, "m^-1"), surface_site_density=(2.72e-9, "mol/cm^2"), - coverage_dependence=True, + thermo_coverage_dependence=True, termination=[], ) # in chemkin, the sites are mostly occupied in about 1e-8 seconds. @@ -1052,13 +1057,21 @@ def test_solve_ch3_thermo_coverage_dependence(self): rxn1.get_surface_rate_coefficient(rxn_system.T.value_si, rxn_system.surface_site_density.value_si), ) - assert isinstance(ch3.thermo.thermo_coverage_dependence, dict) # check to make sure coverage_dependence is still the correct type - for species, parameters in ch3.thermo.thermo_coverage_dependence.items(): - assert isinstance(species, Species) # species should be a Species + assert isinstance(ch3x.thermo.thermo_coverage_dependence, dict) # check to make sure coverage_dependence is still the correct type + for species, parameters in ch3x.thermo.thermo_coverage_dependence.items(): + assert isinstance(species, str) # species should be an ajacency list assert isinstance(parameters, dict) assert parameters["model"] is not None assert parameters["enthalpy-coefficients"] is not None assert parameters["entropy-coefficients"] is not None + + # check thermo_coverage_dependence is on + # and thermo_coeff_matrix and stoi_matrix are correctly created + assert rxn_system.thermo_coverage_dependence is True + assert np.array_equal(rxn_system.stoi_matrix, np.array([[-1, -1, 1]], dtype=float)) + thermo_coeff_matrix = np.array([np.zeros((3,6))]*3) + thermo_coeff_matrix[-1][-1] = [1e5, 2, 3, 1, 5, 3] + assert np.array_equal(rxn_system.thermo_coeff_matrix, thermo_coeff_matrix) # Integrate to get the solution at each time point t = [] @@ -1122,6 +1135,12 @@ def test_solve_ch3_thermo_coverage_dependence(self): ) rxn_system.initialize_model(core_species, core_reactions, edge_species, edge_reactions) + + # check thermo_coverage_dependence is off + # and thermo_coeff_matrix and stoi_matrix are not created + assert rxn_system.thermo_coverage_dependence is False + assert rxn_system.thermo_coeff_matrix is None + assert rxn_system.stoi_matrix is None tlist = np.logspace(-13, -5, 81, dtype=float) From e6fac7b813650930a7153238a92eeffd6bea560b Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 25 Mar 2024 14:12:03 -0400 Subject: [PATCH 20/49] save thermo coverage dependent models as numpy arrays instead a list of rmg quantities --- rmgpy/thermo/nasa.pyx | 5 ++--- rmgpy/thermo/thermodata.pyx | 4 ++-- rmgpy/thermo/wilhoit.pyx | 4 ++-- 3 files changed, 6 insertions(+), 7 deletions(-) diff --git a/rmgpy/thermo/nasa.pyx b/rmgpy/thermo/nasa.pyx index 4014cdfc1be..d749f986cb9 100644 --- a/rmgpy/thermo/nasa.pyx +++ b/rmgpy/thermo/nasa.pyx @@ -34,7 +34,6 @@ cimport numpy as np from libc.math cimport log cimport rmgpy.constants as constants -import rmgpy.quantity as quantity ################################################################################ @@ -314,8 +313,8 @@ cdef class NASA(HeatCapacityModel): for species, parameters in value.items(): # just the polynomial model for now processed_parameters = {'model': parameters['model'], - 'enthalpy-coefficients': [quantity.Dimensionless(p) for p in parameters['enthalpy-coefficients']], - 'entropy-coefficients': [quantity.Dimensionless(p) for p in parameters['entropy-coefficients']], + 'enthalpy-coefficients': np.array([p for p in parameters['enthalpy-coefficients']]), + 'entropy-coefficients': np.array([p for p in parameters['entropy-coefficients']]), } self._thermo_coverage_dependence[species] = processed_parameters diff --git a/rmgpy/thermo/thermodata.pyx b/rmgpy/thermo/thermodata.pyx index e38ae87b51b..9274e088f9a 100644 --- a/rmgpy/thermo/thermodata.pyx +++ b/rmgpy/thermo/thermodata.pyx @@ -128,8 +128,8 @@ cdef class ThermoData(HeatCapacityModel): for species, parameters in value.items(): # just the polynomial model for now processed_parameters = {'model': parameters['model'], - 'enthalpy-coefficients': [quantity.Dimensionless(p) for p in parameters['enthalpy-coefficients']], - 'entropy-coefficients': [quantity.Dimensionless(p) for p in parameters['entropy-coefficients']], + 'enthalpy-coefficients': np.array([p for p in parameters['enthalpy-coefficients']]), + 'entropy-coefficients': np.array([p for p in parameters['entropy-coefficients']]), } self._thermo_coverage_dependence[species] = processed_parameters diff --git a/rmgpy/thermo/wilhoit.pyx b/rmgpy/thermo/wilhoit.pyx index 15299666364..8801511491f 100644 --- a/rmgpy/thermo/wilhoit.pyx +++ b/rmgpy/thermo/wilhoit.pyx @@ -150,8 +150,8 @@ cdef class Wilhoit(HeatCapacityModel): for species, parameters in value.items(): # just the polynomial model for now processed_parameters = {'model': parameters['model'], - 'enthalpy-coefficients': [quantity.Dimensionless(p) for p in parameters['enthalpy-coefficients']], - 'entropy-coefficients': [quantity.Dimensionless(p) for p in parameters['entropy-coefficients']], + 'enthalpy-coefficients': np.array([p for p in parameters['enthalpy-coefficients']]), + 'entropy-coefficients': np.array([p for p in parameters['entropy-coefficients']]), } self._thermo_coverage_dependence[species] = processed_parameters From 0360c5796f42541dabafba9e065d9025165735c3 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 25 Mar 2024 17:47:05 -0400 Subject: [PATCH 21/49] create a subclass of numpy array overriding __repr__ method --- rmgpy/util.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/rmgpy/util.py b/rmgpy/util.py index 15769e02e61..365833379ea 100644 --- a/rmgpy/util.py +++ b/rmgpy/util.py @@ -33,6 +33,7 @@ import shutil import time from functools import wraps +import numpy as np class Subject(object): @@ -238,3 +239,13 @@ def as_list(item, default=None): return default else: return [item] + +class np_list(np.ndarray): + """ + A subclass of numpy.ndarray which rendered as a list when printed. + """ + def __new__(cls, input_array): + obj = np.asarray(input_array).view(cls) + return obj + def __repr__(self): + return str(self.tolist()) From 0c77c57170c9ce735d94ede440501e0778457ebd Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 25 Mar 2024 17:52:59 -0400 Subject: [PATCH 22/49] create thermo coverage dep polynomial model in numpy array class with new __repr__ method --- rmgpy/thermo/nasa.pyx | 5 +++-- rmgpy/thermo/thermodata.pyx | 5 +++-- rmgpy/thermo/wilhoit.pyx | 8 +++++--- 3 files changed, 11 insertions(+), 7 deletions(-) diff --git a/rmgpy/thermo/nasa.pyx b/rmgpy/thermo/nasa.pyx index d749f986cb9..e270400ae5b 100644 --- a/rmgpy/thermo/nasa.pyx +++ b/rmgpy/thermo/nasa.pyx @@ -34,6 +34,7 @@ cimport numpy as np from libc.math cimport log cimport rmgpy.constants as constants +from rmgpy.util import np_list ################################################################################ @@ -313,8 +314,8 @@ cdef class NASA(HeatCapacityModel): for species, parameters in value.items(): # just the polynomial model for now processed_parameters = {'model': parameters['model'], - 'enthalpy-coefficients': np.array([p for p in parameters['enthalpy-coefficients']]), - 'entropy-coefficients': np.array([p for p in parameters['entropy-coefficients']]), + 'enthalpy-coefficients': np_list([p for p in parameters['enthalpy-coefficients']]), + 'entropy-coefficients': np_list([p for p in parameters['entropy-coefficients']]), } self._thermo_coverage_dependence[species] = processed_parameters diff --git a/rmgpy/thermo/thermodata.pyx b/rmgpy/thermo/thermodata.pyx index 9274e088f9a..aed7d18f89a 100644 --- a/rmgpy/thermo/thermodata.pyx +++ b/rmgpy/thermo/thermodata.pyx @@ -35,6 +35,7 @@ cimport numpy as np from libc.math cimport log import rmgpy.quantity as quantity +from rmgpy.util import np_list ################################################################################ @@ -128,8 +129,8 @@ cdef class ThermoData(HeatCapacityModel): for species, parameters in value.items(): # just the polynomial model for now processed_parameters = {'model': parameters['model'], - 'enthalpy-coefficients': np.array([p for p in parameters['enthalpy-coefficients']]), - 'entropy-coefficients': np.array([p for p in parameters['entropy-coefficients']]), + 'enthalpy-coefficients': np_list([p for p in parameters['enthalpy-coefficients']]), + 'entropy-coefficients': np_list([p for p in parameters['entropy-coefficients']]), } self._thermo_coverage_dependence[species] = processed_parameters diff --git a/rmgpy/thermo/wilhoit.pyx b/rmgpy/thermo/wilhoit.pyx index 8801511491f..163b39e1876 100644 --- a/rmgpy/thermo/wilhoit.pyx +++ b/rmgpy/thermo/wilhoit.pyx @@ -35,6 +35,7 @@ from libc.math cimport sqrt, log cimport rmgpy.constants as constants import rmgpy.quantity as quantity +from rmgpy.util import np_list # Prior to numpy 1.14, `numpy.linalg.lstsq` does not accept None as a value RCOND = -1 if int(np.__version__.split('.')[1]) < 14 else None @@ -150,8 +151,8 @@ cdef class Wilhoit(HeatCapacityModel): for species, parameters in value.items(): # just the polynomial model for now processed_parameters = {'model': parameters['model'], - 'enthalpy-coefficients': np.array([p for p in parameters['enthalpy-coefficients']]), - 'entropy-coefficients': np.array([p for p in parameters['entropy-coefficients']]), + 'enthalpy-coefficients': np_list([p for p in parameters['enthalpy-coefficients']]), + 'entropy-coefficients': np_list([p for p in parameters['entropy-coefficients']]), } self._thermo_coverage_dependence[species] = processed_parameters @@ -213,7 +214,8 @@ cdef class Wilhoit(HeatCapacityModel): self.Cp0, self.CpInf, self.a0, self.a1, self.a2, self.a3, self.H0, self.S0, self.B, - Tmin=self.Tmin, Tmax=self.Tmax, comment=self.comment, + Tmin=self.Tmin, Tmax=self.Tmax, thermo_coverage_dependence=self.thermo_coverage_dependence, + comment=self.comment, ) @cython.boundscheck(False) From 4f09c50fbb3b817617ef5b101213a544c522e960 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 25 Mar 2024 17:53:56 -0400 Subject: [PATCH 23/49] use adjacency list as the key of thermo cov dep models --- test/rmgpy/thermo/nasaTest.py | 9 ++++++--- test/rmgpy/thermo/thermodataTest.py | 2 +- test/rmgpy/thermo/wilhoitTest.py | 19 +++++++------------ 3 files changed, 14 insertions(+), 16 deletions(-) diff --git a/test/rmgpy/thermo/nasaTest.py b/test/rmgpy/thermo/nasaTest.py index 8945bef5f42..f96df94cfe0 100644 --- a/test/rmgpy/thermo/nasaTest.py +++ b/test/rmgpy/thermo/nasaTest.py @@ -39,7 +39,6 @@ import rmgpy.constants as constants from rmgpy.quantity import ScalarQuantity from rmgpy.thermo.nasa import NASA, NASAPolynomial -from rmgpy.quantity import Dimensionless class TestNASA: @@ -70,7 +69,7 @@ def setup_class(self): self.Tmax = 3000.0 self.Tint = 650.73 self.E0 = -782292.0 # J/mol. - self.thermo_coverage_dependence = {'OX':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}} + self.thermo_coverage_dependence = {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}} self.comment = "C2H6" self.nasa = NASA( polynomials=[ @@ -377,7 +376,11 @@ def test_nasa_as_dict_full(self): assert nasa_dict["E0"]["value"] == self.E0 assert nasa_dict["Tmin"]["value"] == self.Tmin assert nasa_dict["Tmax"]["value"] == self.Tmax - assert repr(nasa_dict["thermo_coverage_dependence"]) == "{'OX': {'model': 'polynomial', 'enthalpy-coefficients': [{'class': 'ScalarQuantity', 'value': 1.0}, {'class': 'ScalarQuantity', 'value': 2.0}, {'class': 'ScalarQuantity', 'value': 3.0}], 'entropy-coefficients': [{'class': 'ScalarQuantity', 'value': 1.0}, {'class': 'ScalarQuantity', 'value': 2.0}, {'class': 'ScalarQuantity', 'value': 3.0}]}}" + assert nasa_dict["thermo_coverage_dependence"].keys() == self.thermo_coverage_dependence.keys() + sp_name = list(self.thermo_coverage_dependence.keys())[0] + assert nasa_dict['thermo_coverage_dependence'][sp_name]['model'] == self.thermo_coverage_dependence[sp_name]['model'] + assert nasa_dict['thermo_coverage_dependence'][sp_name]['enthalpy-coefficients']['object'] == self.thermo_coverage_dependence[sp_name]['enthalpy-coefficients'] + assert nasa_dict['thermo_coverage_dependence'][sp_name]['entropy-coefficients']['object'] == self.thermo_coverage_dependence[sp_name]['entropy-coefficients'] assert nasa_dict["comment"] == self.comment assert tuple(nasa_dict["polynomials"]["polynomial1"]["coeffs"]["object"]) == tuple(self.coeffs_low) assert tuple(nasa_dict["polynomials"]["polynomial2"]["coeffs"]["object"]) == tuple(self.coeffs_high) diff --git a/test/rmgpy/thermo/thermodataTest.py b/test/rmgpy/thermo/thermodataTest.py index cca336e1332..3c6b39e7e03 100644 --- a/test/rmgpy/thermo/thermodataTest.py +++ b/test/rmgpy/thermo/thermodataTest.py @@ -53,7 +53,7 @@ def setup_class(self): self.Tmin = 100.0 self.Tmax = 3000.0 self.E0 = -782292.0 - self.thermo_coverage_dependence = {'OX':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}} + self.thermo_coverage_dependence = {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}} self.comment = "C2H6" self.thermodata = ThermoData( Tdata=(self.Tdata, "K"), diff --git a/test/rmgpy/thermo/wilhoitTest.py b/test/rmgpy/thermo/wilhoitTest.py index 15d76d34fa1..6a350d78336 100644 --- a/test/rmgpy/thermo/wilhoitTest.py +++ b/test/rmgpy/thermo/wilhoitTest.py @@ -58,7 +58,7 @@ def setup_class(self): self.Tmin = 300.0 self.Tmax = 3000.0 self.comment = "C2H6" - self.thermo_coverage_dependence = {'OX':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}} + self.thermo_coverage_dependence = {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}} self.wilhoit = Wilhoit( Cp0=(self.Cp0 * constants.R, "J/(mol*K)"), CpInf=(self.CpInf * constants.R, "J/(mol*K)"), @@ -458,16 +458,11 @@ def test_wilhoit_as_dict(self): "value": 178.76114800000002, }, "class": "Wilhoit", - 'thermo_coverage_dependence': {'OX': { - 'model': 'polynomial', - 'enthalpy-coefficients': [{'class': 'ScalarQuantity', 'value': 1.0}, - {'class': 'ScalarQuantity', 'value': 2.0}, - {'class': 'ScalarQuantity', 'value': 3.0} - ], - 'entropy-coefficients': [{'class': 'ScalarQuantity', 'value': 1.0}, - {'class': 'ScalarQuantity', 'value': 2.0}, - {'class': 'ScalarQuantity', 'value': 3.0} - ]}} + 'thermo_coverage_dependence': {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}': { + 'model': 'polynomial', + 'enthalpy-coefficients': {'class': 'np_array', 'object': [1, 2, 3]}, + 'entropy-coefficients': {'class': 'np_array', 'object': [1, 2, 3]}} + } } def test_make_wilhoit(self): @@ -476,6 +471,6 @@ def test_make_wilhoit(self): """ wilhoit_dict = self.wilhoit.as_dict() new_wilhoit = Wilhoit.__new__(Wilhoit) - class_dictionary = {"ScalarQuantity": ScalarQuantity, "Wilhoit": Wilhoit} + class_dictionary = {"ScalarQuantity": ScalarQuantity, "Wilhoit": Wilhoit, "np_array": np.array} new_wilhoit.make_object(wilhoit_dict, class_dictionary) From 31de2b18c29c206eda69201bab9243f6b923a50d Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 25 Mar 2024 18:01:32 -0400 Subject: [PATCH 24/49] adjust the method to read data from numpy arrays instead list --- rmgpy/solver/surface.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index becdd844f46..1c602734063 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -215,8 +215,8 @@ cdef class SurfaceReactor(ReactionSystem): for species in self.species_index.keys(): if species.is_isomorphic(molecule, strict=False): species_index = self.species_index[species] - thermo_polynomials = parameters['enthalpy-coefficients'] + parameters['entropy-coefficients'] - self.thermo_coeff_matrix[sp_index, species_index] = [x.value_si for x in thermo_polynomials] + thermo_polynomials = np.concatenate((parameters['enthalpy-coefficients'], parameters['entropy-coefficients']), 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: From 43dea6bc027db0900776b995c756593db9d8e803 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Thu, 28 Mar 2024 14:16:08 -0400 Subject: [PATCH 25/49] Write thermo coverage dependent data to Cantera yaml file --- rmgpy/rmg/main.py | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index 66a827e8890..a7a7ac63cf0 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -49,6 +49,7 @@ import yaml from cantera import ck2yaml from scipy.optimize import brute +import cantera as ct import rmgpy.util as util from rmgpy import settings @@ -1229,6 +1230,42 @@ 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: + # add thermo coverage dependence to Cantera files + chem_yaml_path = os.path.join(self.output_directory, "cantera", "chem.yaml") + gas = ct.Solution(chem_yaml_path, "gas") + surf = ct.Interface(chem_yaml_path, "surface1", [gas]) + with open(chem_yaml_path, 'r') as f: + content = yaml.load(f, Loader=yaml.FullLoader) + + content['phases'][1]['reference-state-coverage'] = 0.11 + content['phases'][1]['thermo'] = 'coverage-dependent-surface' + + for s in self.reaction_model.core.species: + if s.contains_surface_site() and s.thermo.thermo_coverage_dependence: + for dep_sp, parameters in s.thermo.thermo_coverage_dependence.items(): + mol = Molecule().from_adjacency_list(dep_sp) + for sp in self.reaction_model.core.species: + if sp.is_isomorphic(mol, strict=False): + try: + parameters['units'] = {'energy':'J', 'quantity':'mol'} + content["species"][surf.species_index(sp.label)]['coverage-dependencies'][sp.label] = parameters + except KeyError: + content["species"][surf.species_index(sp.label)]['coverage-dependencies'] = {sp.label: parameters} + + annotated_yaml_path = os.path.join(self.output_directory, "cantera", "chem_annotated.yaml") + with open(annotated_yaml_path, 'r') as f: + annotated_content = yaml.load(f, Loader=yaml.FullLoader) + + annotated_content['phases'] = content['phases'] + annotated_content['species'] = content['species'] + + with open(chem_yaml_path, 'w') as output_f: + yaml.dump(content, output_f, sort_keys=False) + + with open(annotated_yaml_path, 'w') as output_f: + yaml.dump(annotated_content, output_f, sort_keys=False) + 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")) From ed0474cdf5249324d2f9656c634a89eba9501e72 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Fri, 29 Mar 2024 14:29:47 -0400 Subject: [PATCH 26/49] Use chemkin strings instead of labels when writing to yaml files --- rmgpy/rmg/main.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index a7a7ac63cf0..542d015cdfe 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -1249,9 +1249,9 @@ def execute(self, initialize=True, **kwargs): if sp.is_isomorphic(mol, strict=False): try: parameters['units'] = {'energy':'J', 'quantity':'mol'} - content["species"][surf.species_index(sp.label)]['coverage-dependencies'][sp.label] = parameters + content["species"][surf.species_index(sp.to_chemkin())]['coverage-dependencies'][sp.to_chemkin()] = parameters except KeyError: - content["species"][surf.species_index(sp.label)]['coverage-dependencies'] = {sp.label: parameters} + content["species"][surf.species_index(sp.to_chemkin())]['coverage-dependencies'] = {sp.to_chemkin(): parameters} annotated_yaml_path = os.path.join(self.output_directory, "cantera", "chem_annotated.yaml") with open(annotated_yaml_path, 'r') as f: From 51e3821185a1616821d4a1d3f3d40cbb579f34d6 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Fri, 29 Mar 2024 21:06:43 -0400 Subject: [PATCH 27/49] Write the thermo coverage dependent data to yaml files in right format --- rmgpy/rmg/main.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index 542d015cdfe..e1826b95da9 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -1247,11 +1247,13 @@ def execute(self, initialize=True, **kwargs): mol = Molecule().from_adjacency_list(dep_sp) for sp in self.reaction_model.core.species: if sp.is_isomorphic(mol, strict=False): + parameters['units'] = {'energy':'J', 'quantity':'mol'} + parameters['enthalpy-coefficients'] = [float(value) for value in parameters['enthalpy-coefficients']] + parameters['entropy-coefficients'] = [float(value) for value in parameters['entropy-coefficients']] try: - parameters['units'] = {'energy':'J', 'quantity':'mol'} - content["species"][surf.species_index(sp.to_chemkin())]['coverage-dependencies'][sp.to_chemkin()] = parameters + content["species"][gas.n_species+surf.species_index(sp.to_chemkin())]['coverage-dependencies'][sp.to_chemkin()] = parameters except KeyError: - content["species"][surf.species_index(sp.to_chemkin())]['coverage-dependencies'] = {sp.to_chemkin(): parameters} + content["species"][gas.n_species+surf.species_index(sp.to_chemkin())]['coverage-dependencies'] = {sp.to_chemkin(): parameters} annotated_yaml_path = os.path.join(self.output_directory, "cantera", "chem_annotated.yaml") with open(annotated_yaml_path, 'r') as f: From 79098edcd1a6fc90707242bdd5a29140ad036079 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Tue, 2 Apr 2024 12:50:51 -0400 Subject: [PATCH 28/49] Add thermo_coverag_dependence as an instance for class RMG --- rmgpy/rmg/main.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index e1826b95da9..c85224f1309 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -194,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 @@ -511,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 @@ -769,7 +770,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: From fae856c30e713fc34f6db7b37dbd276a68b225fc Mon Sep 17 00:00:00 2001 From: 12Chao Date: Wed, 17 Jul 2024 20:40:46 -0400 Subject: [PATCH 29/49] fix the bug that the coverage kinetics cannot be written to chemkin --- rmgpy/chemkin.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index 1988419057d..8b263fd0186 100644 --- a/rmgpy/chemkin.pyx +++ b/rmgpy/chemkin.pyx @@ -1867,7 +1867,7 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals for species, cov_params in kinetics.coverage_dependence.items(): label = get_species_identifier(species) string += f' COV / {label:<41} ' - string += f"{cov_params['a']:<9.3g} {cov_params['m']:<9.3g} {cov_params['E'].value_si/4184.:<9.3f} /\n" + string += f"{cov_params['a'].value_si:<9.3g} {cov_params['m'].value_si:<9.3g} {cov_params['E'].value_si/4184.:<9.3f} /\n" if isinstance(kinetics, (_kinetics.ThirdBody, _kinetics.Lindemann, _kinetics.Troe)): # Write collider efficiencies From b35ee5b8b7f5d3bc937f851538d356c21fbc4018 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Sat, 19 Oct 2024 10:18:17 -0400 Subject: [PATCH 30/49] write data as RMG object instead of floats to keep the units consistent --- rmgpy/rmg/main.py | 4 ++-- rmgpy/solver/surface.pyx | 4 +++- rmgpy/thermo/nasa.pyx | 5 +++-- rmgpy/thermo/thermodata.pyx | 4 ++-- rmgpy/thermo/wilhoit.pyx | 4 ++-- test/rmgpy/thermo/nasaTest.py | 15 ++++++++------- 6 files changed, 20 insertions(+), 16 deletions(-) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index c85224f1309..68940fbc292 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -1249,8 +1249,8 @@ def execute(self, initialize=True, **kwargs): for sp in self.reaction_model.core.species: if sp.is_isomorphic(mol, strict=False): parameters['units'] = {'energy':'J', 'quantity':'mol'} - parameters['enthalpy-coefficients'] = [float(value) for value in parameters['enthalpy-coefficients']] - parameters['entropy-coefficients'] = [float(value) for value in parameters['entropy-coefficients']] + parameters['enthalpy-coefficients'] = [value.value_si for value in parameters['enthalpy-coefficients']] + parameters['entropy-coefficients'] = [value.value_si for value in parameters['entropy-coefficients']] try: content["species"][gas.n_species+surf.species_index(sp.to_chemkin())]['coverage-dependencies'][sp.to_chemkin()] = parameters except KeyError: diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index 1c602734063..2fad5164f16 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -215,7 +215,9 @@ cdef class SurfaceReactor(ReactionSystem): for species in self.species_index.keys(): if species.is_isomorphic(molecule, strict=False): species_index = self.species_index[species] - thermo_polynomials = np.concatenate((parameters['enthalpy-coefficients'], parameters['entropy-coefficients']), axis=0) + 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 diff --git a/rmgpy/thermo/nasa.pyx b/rmgpy/thermo/nasa.pyx index e270400ae5b..0b04438654a 100644 --- a/rmgpy/thermo/nasa.pyx +++ b/rmgpy/thermo/nasa.pyx @@ -35,6 +35,7 @@ from libc.math cimport log cimport rmgpy.constants as constants from rmgpy.util import np_list +import rmgpy.quantity as quantity ################################################################################ @@ -314,8 +315,8 @@ cdef class NASA(HeatCapacityModel): for species, parameters in value.items(): # just the polynomial model for now processed_parameters = {'model': parameters['model'], - 'enthalpy-coefficients': np_list([p for p in parameters['enthalpy-coefficients']]), - 'entropy-coefficients': np_list([p for p in parameters['entropy-coefficients']]), + 'enthalpy-coefficients': np_list([quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']]), + 'entropy-coefficients': np_list([quantity.Entropy(p) for p in parameters['entropy-coefficients']]), } self._thermo_coverage_dependence[species] = processed_parameters diff --git a/rmgpy/thermo/thermodata.pyx b/rmgpy/thermo/thermodata.pyx index aed7d18f89a..6797fa99e4d 100644 --- a/rmgpy/thermo/thermodata.pyx +++ b/rmgpy/thermo/thermodata.pyx @@ -129,8 +129,8 @@ cdef class ThermoData(HeatCapacityModel): for species, parameters in value.items(): # just the polynomial model for now processed_parameters = {'model': parameters['model'], - 'enthalpy-coefficients': np_list([p for p in parameters['enthalpy-coefficients']]), - 'entropy-coefficients': np_list([p for p in parameters['entropy-coefficients']]), + 'enthalpy-coefficients': np_list([quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']]), + 'entropy-coefficients': np_list([quantity.Entropy(p) for p in parameters['entropy-coefficients']]), } self._thermo_coverage_dependence[species] = processed_parameters diff --git a/rmgpy/thermo/wilhoit.pyx b/rmgpy/thermo/wilhoit.pyx index 163b39e1876..3342f64c48f 100644 --- a/rmgpy/thermo/wilhoit.pyx +++ b/rmgpy/thermo/wilhoit.pyx @@ -151,8 +151,8 @@ cdef class Wilhoit(HeatCapacityModel): for species, parameters in value.items(): # just the polynomial model for now processed_parameters = {'model': parameters['model'], - 'enthalpy-coefficients': np_list([p for p in parameters['enthalpy-coefficients']]), - 'entropy-coefficients': np_list([p for p in parameters['entropy-coefficients']]), + 'enthalpy-coefficients': np_list([quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']]), + 'entropy-coefficients': np_list([quantity.Entropy(p) for p in parameters['entropy-coefficients']]), } self._thermo_coverage_dependence[species] = processed_parameters diff --git a/test/rmgpy/thermo/nasaTest.py b/test/rmgpy/thermo/nasaTest.py index f96df94cfe0..23255082f9a 100644 --- a/test/rmgpy/thermo/nasaTest.py +++ b/test/rmgpy/thermo/nasaTest.py @@ -31,8 +31,7 @@ This script contains unit tests of the :mod:`rmgpy.thermo.nasa` module. """ -import os.path - +import os.path, ast import numpy as np @@ -69,7 +68,7 @@ def setup_class(self): self.Tmax = 3000.0 self.Tint = 650.73 self.E0 = -782292.0 # J/mol. - self.thermo_coverage_dependence = {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}} + self.thermo_coverage_dependence = {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}':{'model':'polynomial', 'enthalpy-coefficients':[(1,'J/mol'),(2,'J/mol'),(3,'J/mol')], "entropy-coefficients":[(1,'J/(mol*K)'),(2,'J/(mol*K)'),(3,'J/(mol*K)')]}} self.comment = "C2H6" self.nasa = NASA( polynomials=[ @@ -143,7 +142,7 @@ def test_thermo_coverage_dependence(self): """ Test that the thermo_coverage_dependence property was properly set. """ - assert repr(self.nasa.thermo_coverage_dependence) == repr(self.thermo_coverage_dependence) + assert ast.literal_eval(repr(self.nasa.thermo_coverage_dependence)) == ast.literal_eval(repr(self.thermo_coverage_dependence)) def test_is_temperature_valid(self): """ @@ -271,7 +270,7 @@ def test_pickle(self): assert self.nasa.Tmax.units == nasa.Tmax.units assert self.nasa.E0.value == nasa.E0.value assert self.nasa.E0.units == nasa.E0.units - assert repr(self.nasa.thermo_coverage_dependence) == repr(nasa.thermo_coverage_dependence) + assert ast.literal_eval(repr(self.nasa.thermo_coverage_dependence)) == ast.literal_eval(repr(nasa.thermo_coverage_dependence)) assert self.nasa.comment == nasa.comment def test_repr(self): @@ -305,7 +304,7 @@ def test_repr(self): assert self.nasa.Tmax.units == nasa.Tmax.units assert self.nasa.E0.value == nasa.E0.value assert self.nasa.E0.units == nasa.E0.units - assert repr(self.nasa.thermo_coverage_dependence) == repr(nasa.thermo_coverage_dependence) + assert ast.literal_eval(repr(self.nasa.thermo_coverage_dependence)) == ast.literal_eval(repr(nasa.thermo_coverage_dependence)) assert self.nasa.comment == nasa.comment def test_to_cantera(self): @@ -379,7 +378,9 @@ def test_nasa_as_dict_full(self): assert nasa_dict["thermo_coverage_dependence"].keys() == self.thermo_coverage_dependence.keys() sp_name = list(self.thermo_coverage_dependence.keys())[0] assert nasa_dict['thermo_coverage_dependence'][sp_name]['model'] == self.thermo_coverage_dependence[sp_name]['model'] - assert nasa_dict['thermo_coverage_dependence'][sp_name]['enthalpy-coefficients']['object'] == self.thermo_coverage_dependence[sp_name]['enthalpy-coefficients'] + enthalpy_list = nasa_dict['thermo_coverage_dependence'][sp_name]['enthalpy-coefficients']['object'] + # return [(str(coeff.value), str(coeff.units))for coeff in enthalpy_list], self.thermo_coverage_dependence[sp_name]['enthalpy-coefficients'] + assert [(int(coeff.value), str(coeff.units))for coeff in enthalpy_list] == self.thermo_coverage_dependence[sp_name]['enthalpy-coefficients'] assert nasa_dict['thermo_coverage_dependence'][sp_name]['entropy-coefficients']['object'] == self.thermo_coverage_dependence[sp_name]['entropy-coefficients'] assert nasa_dict["comment"] == self.comment assert tuple(nasa_dict["polynomials"]["polynomial1"]["coeffs"]["object"]) == tuple(self.coeffs_low) From 2ac9c17c95ca529fa20c26260e5e180ce4838578 Mon Sep 17 00:00:00 2001 From: Bjarne Kreitz Date: Fri, 10 Oct 2025 13:26:08 -0400 Subject: [PATCH 31/49] add thermo_coverage_dependence to initialize --- rmgpy/rmg/main.py | 1 + 1 file changed, 1 insertion(+) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index 68940fbc292..98a027fff12 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -525,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 From a579b3ab6ec2108eedb00d5c0eb7d5761337bb1d Mon Sep 17 00:00:00 2001 From: Bjarne Kreitz Date: Mon, 13 Oct 2025 14:58:34 -0400 Subject: [PATCH 32/49] fix type in wilhoit.pyx --- rmgpy/thermo/wilhoit.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/thermo/wilhoit.pyx b/rmgpy/thermo/wilhoit.pyx index 3342f64c48f..bcbf16d4abb 100644 --- a/rmgpy/thermo/wilhoit.pyx +++ b/rmgpy/thermo/wilhoit.pyx @@ -59,7 +59,7 @@ cdef class Wilhoit(HeatCapacityModel): `B` The Wilhoit scaled temperature coefficient in K `Tmin` The minimum temperature in K at which the model is valid, or zero if unknown or undefined `Tmax` The maximum temperature in K at which the model is valid, or zero if unknown or undefined - `thermo_covreage_dependence` The coverage dependence of the thermo + `thermo_coverage_dependence` The coverage dependence of the thermo `comment` Information about the model (e.g. its source) ============================= ========================================================================================= """ From b504d9ecbee229776d1c426335583d639ee3d487 Mon Sep 17 00:00:00 2001 From: Bjarne Kreitz Date: Mon, 13 Oct 2025 17:22:54 -0400 Subject: [PATCH 33/49] fix coefficients check in nasaTest.py --- test/rmgpy/thermo/nasaTest.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/rmgpy/thermo/nasaTest.py b/test/rmgpy/thermo/nasaTest.py index 23255082f9a..138ea59d3f7 100644 --- a/test/rmgpy/thermo/nasaTest.py +++ b/test/rmgpy/thermo/nasaTest.py @@ -378,10 +378,10 @@ def test_nasa_as_dict_full(self): assert nasa_dict["thermo_coverage_dependence"].keys() == self.thermo_coverage_dependence.keys() sp_name = list(self.thermo_coverage_dependence.keys())[0] assert nasa_dict['thermo_coverage_dependence'][sp_name]['model'] == self.thermo_coverage_dependence[sp_name]['model'] - enthalpy_list = nasa_dict['thermo_coverage_dependence'][sp_name]['enthalpy-coefficients']['object'] - # return [(str(coeff.value), str(coeff.units))for coeff in enthalpy_list], self.thermo_coverage_dependence[sp_name]['enthalpy-coefficients'] - assert [(int(coeff.value), str(coeff.units))for coeff in enthalpy_list] == self.thermo_coverage_dependence[sp_name]['enthalpy-coefficients'] - assert nasa_dict['thermo_coverage_dependence'][sp_name]['entropy-coefficients']['object'] == self.thermo_coverage_dependence[sp_name]['entropy-coefficients'] + enthalpy_list = nasa_dict['thermo_coverage_dependence'][sp_name]['enthalpy-coefficients']['object'] + assert [(int(coeff.value), str(coeff.units)) for coeff in enthalpy_list] == self.thermo_coverage_dependence[sp_name]['enthalpy-coefficients'] + entropy_list = nasa_dict['thermo_coverage_dependence'][sp_name]['entropy-coefficients']['object'] + assert [(int(coeff.value), str(coeff.units)) for coeff in entropy_list] == self.thermo_coverage_dependence[sp_name]['entropy-coefficients'] assert nasa_dict["comment"] == self.comment assert tuple(nasa_dict["polynomials"]["polynomial1"]["coeffs"]["object"]) == tuple(self.coeffs_low) assert tuple(nasa_dict["polynomials"]["polynomial2"]["coeffs"]["object"]) == tuple(self.coeffs_high) From 91fbefb8f883931a9cc95e19ecfc3d5d8ef464ce Mon Sep 17 00:00:00 2001 From: Bjarne Kreitz Date: Mon, 13 Oct 2025 20:57:29 -0400 Subject: [PATCH 34/49] fix test_thermo_coverage_dependence --- test/rmgpy/thermo/wilhoitTest.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/rmgpy/thermo/wilhoitTest.py b/test/rmgpy/thermo/wilhoitTest.py index 6a350d78336..b3166796025 100644 --- a/test/rmgpy/thermo/wilhoitTest.py +++ b/test/rmgpy/thermo/wilhoitTest.py @@ -39,7 +39,7 @@ import rmgpy.constants as constants from rmgpy.quantity import ScalarQuantity from rmgpy.thermo.wilhoit import Wilhoit - +import rmgpy.quantity as quantity class TestWilhoit: """ @@ -58,7 +58,7 @@ def setup_class(self): self.Tmin = 300.0 self.Tmax = 3000.0 self.comment = "C2H6" - self.thermo_coverage_dependence = {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}} + self.thermo_coverage_dependence = {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}':{'model':'polynomial', 'enthalpy-coefficients':[quantity.Enthalpy(1,'J/mol'),quantity.Enthalpy(2,'J/mol'),quantity.Enthalpy(3,'J/mol')], "entropy-coefficients":[quantity.Entropy(1,'J/(mol*K)'),quantity.Entropy(2,'J/(mol*K)'),quantity.Entropy(3,'J/(mol*K)')]}} self.wilhoit = Wilhoit( Cp0=(self.Cp0 * constants.R, "J/(mol*K)"), CpInf=(self.CpInf * constants.R, "J/(mol*K)"), @@ -460,8 +460,8 @@ def test_wilhoit_as_dict(self): "class": "Wilhoit", 'thermo_coverage_dependence': {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}': { 'model': 'polynomial', - 'enthalpy-coefficients': {'class': 'np_array', 'object': [1, 2, 3]}, - 'entropy-coefficients': {'class': 'np_array', 'object': [1, 2, 3]}} + 'enthalpy-coefficients': {'class': 'np_array', 'object': [(1,'J/mol'), (2,'J/mol'), (3,'J/mol')]}, + 'entropy-coefficients': {'class': 'np_array', 'object': [(1,'J/(mol*K)'),(2,'J/(mol*K)'),(3,'J/(mol*K)')]}} } } From e5ecc8d4c28d1a253e61e01fd8bcc3f665729622 Mon Sep 17 00:00:00 2001 From: Bjarne Kreitz Date: Mon, 13 Oct 2025 20:59:25 -0400 Subject: [PATCH 35/49] fix thermodataTest.py --- test/rmgpy/thermo/thermodataTest.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/rmgpy/thermo/thermodataTest.py b/test/rmgpy/thermo/thermodataTest.py index 3c6b39e7e03..3a0abed4b92 100644 --- a/test/rmgpy/thermo/thermodataTest.py +++ b/test/rmgpy/thermo/thermodataTest.py @@ -36,7 +36,7 @@ import rmgpy.constants as constants from rmgpy.thermo.thermodata import ThermoData - +import rmgpy.quantity as quantity class TestThermoData: """ @@ -53,7 +53,7 @@ def setup_class(self): self.Tmin = 100.0 self.Tmax = 3000.0 self.E0 = -782292.0 - self.thermo_coverage_dependence = {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}} + self.thermo_coverage_dependence = {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}':{'model':'polynomial', 'enthalpy-coefficients':[quantity.Enthalpy(1,'J/mol'),quantity.Enthalpy(2,'J/mol'),quantity.Enthalpy(3,'J/mol')], "entropy-coefficients":[quantity.Entropy(1,'J/(mol*K)'),quantity.Entropy(2,'J/(mol*K)'),quantity.Entropy(3,'J/(mol*K)')]}} self.comment = "C2H6" self.thermodata = ThermoData( Tdata=(self.Tdata, "K"), From a413a2854aef9b6454be78e2b8c6f5a8ea361675 Mon Sep 17 00:00:00 2001 From: Bjarne Kreitz Date: Mon, 13 Oct 2025 21:14:42 -0400 Subject: [PATCH 36/49] fix convertTest.py --- test/rmgpy/thermo/convertTest.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/test/rmgpy/thermo/convertTest.py b/test/rmgpy/thermo/convertTest.py index ef2f7265e31..a60e0ddf606 100644 --- a/test/rmgpy/thermo/convertTest.py +++ b/test/rmgpy/thermo/convertTest.py @@ -38,7 +38,6 @@ import rmgpy.constants as constants from rmgpy.thermo import Wilhoit, NASA, NASAPolynomial, ThermoData - class TestConverter: """ Contains unit tests of the thermodynamics model conversion functions. @@ -57,7 +56,7 @@ def setup_class(self): S0=(-118.46 * constants.R, "J/(mol*K)"), Tmin=(10, "K"), Tmax=(3000, "K"), - thermo_coverage_dependence = {'OX':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}}, + thermo_coverage_dependence = {'OX':{'model':'polynomial', 'enthalpy-coefficients':[(1,'J/mol'),(2,'J/mol'),(3,'J/mol')], "entropy-coefficients":[(1,'J/(mol*K)'),(2,'J/(mol*K)'),(3,'J/(mol*K)')]}}, comment="C2H6", ) self.nasa = NASA( @@ -94,7 +93,7 @@ def setup_class(self): E0=(-93.6077, "kJ/mol"), Cp0=(4.0 * constants.R, "J/(mol*K)"), CpInf=(21.5 * constants.R, "J/(mol*K)"), - thermo_coverage_dependence = {'OX':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}}, + thermo_coverage_dependence = {'OX':{'model':'polynomial', 'enthalpy-coefficients':[(1,'J/mol'),(2,'J/mol'),(3,'J/mol')], "entropy-coefficients":[(1,'J/(mol*K)'),(2,'J/(mol*K)'),(3,'J/(mol*K)')]}}, comment="C2H6", ) self.thermodata = ThermoData( @@ -110,7 +109,7 @@ def setup_class(self): Tmin=(10, "K"), Tmax=(3000, "K"), E0=(-93.6077, "kJ/mol"), - thermo_coverage_dependence = {'OX':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}}, + thermo_coverage_dependence = {'OX':{'model':'polynomial', 'enthalpy-coefficients':[(1,'J/mol'),(2,'J/mol'),(3,'J/mol')], "entropy-coefficients":[(1,'J/(mol*K)'),(2,'J/(mol*K)'),(3,'J/(mol*K)')]}}, comment="C2H6", ) From 6964658a0fae200a44383b82ef284a50473feb6b Mon Sep 17 00:00:00 2001 From: Bjarne Kreitz Date: Tue, 14 Oct 2025 10:08:37 -0400 Subject: [PATCH 37/49] fix wilhoit test_thermo_coverage_dependence and test_wilhoit_as_dict --- test/rmgpy/thermo/wilhoitTest.py | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/test/rmgpy/thermo/wilhoitTest.py b/test/rmgpy/thermo/wilhoitTest.py index b3166796025..92fe4737480 100644 --- a/test/rmgpy/thermo/wilhoitTest.py +++ b/test/rmgpy/thermo/wilhoitTest.py @@ -31,7 +31,7 @@ This script contains unit tests of the :mod:`rmgpy.thermo.wilhoit` module. """ -import os.path +import os.path, ast import numpy as np @@ -58,7 +58,7 @@ def setup_class(self): self.Tmin = 300.0 self.Tmax = 3000.0 self.comment = "C2H6" - self.thermo_coverage_dependence = {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}':{'model':'polynomial', 'enthalpy-coefficients':[quantity.Enthalpy(1,'J/mol'),quantity.Enthalpy(2,'J/mol'),quantity.Enthalpy(3,'J/mol')], "entropy-coefficients":[quantity.Entropy(1,'J/(mol*K)'),quantity.Entropy(2,'J/(mol*K)'),quantity.Entropy(3,'J/(mol*K)')]}} + self.thermo_coverage_dependence = {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}':{'model':'polynomial', 'enthalpy-coefficients':[(1,'J/mol'),(2,'J/mol'),(3,'J/mol')], "entropy-coefficients":[(1,'J/(mol*K)'),(2,'J/(mol*K)'),(3,'J/(mol*K)')]}} self.wilhoit = Wilhoit( Cp0=(self.Cp0 * constants.R, "J/(mol*K)"), CpInf=(self.CpInf * constants.R, "J/(mol*K)"), @@ -165,7 +165,7 @@ def test_thermo_coverage_dependence(self): """ Test that the Wilhoit thermo_coverage_dependence property was properly set. """ - assert repr(self.wilhoit.thermo_coverage_dependence) == repr(self.thermo_coverage_dependence) + assert ast.literal_eval(repr(self.wilhoit.thermo_coverage_dependence)) == ast.literal_eval(repr(self.thermo_coverage_dependence)) def test_is_temperature_valid(self): """ @@ -428,6 +428,16 @@ def test_wilhoit_as_dict(self): Test that a Wilhoit object can be converted to a dictionary representation properly """ wilhoit_dict = self.wilhoit.as_dict() + + assert wilhoit_dict["thermo_coverage_dependence"].keys() == self.thermo_coverage_dependence.keys() + sp_name = list(self.thermo_coverage_dependence.keys())[0] + assert wilhoit_dict['thermo_coverage_dependence'][sp_name]['model'] == self.thermo_coverage_dependence[sp_name]['model'] + enthalpy_list = wilhoit_dict['thermo_coverage_dependence'][sp_name]['enthalpy-coefficients']['object'] + assert [(int(coeff.value), str(coeff.units)) for coeff in enthalpy_list] == self.thermo_coverage_dependence[sp_name]['enthalpy-coefficients'] + entropy_list = wilhoit_dict['thermo_coverage_dependence'][sp_name]['entropy-coefficients']['object'] + assert [(int(coeff.value), str(coeff.units)) for coeff in entropy_list] == self.thermo_coverage_dependence[sp_name]['entropy-coefficients'] + + wilhoit_dict = {k: wilhoit_dict[k] for k in wilhoit_dict.keys() - {'thermo_coverage_dependence'}} assert wilhoit_dict == { "comment": "C2H6", "B": {"units": "K", "class": "ScalarQuantity", "value": 1068.68}, @@ -458,11 +468,11 @@ def test_wilhoit_as_dict(self): "value": 178.76114800000002, }, "class": "Wilhoit", - 'thermo_coverage_dependence': {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}': { - 'model': 'polynomial', - 'enthalpy-coefficients': {'class': 'np_array', 'object': [(1,'J/mol'), (2,'J/mol'), (3,'J/mol')]}, - 'entropy-coefficients': {'class': 'np_array', 'object': [(1,'J/(mol*K)'),(2,'J/(mol*K)'),(3,'J/(mol*K)')]}} - } + # 'thermo_coverage_dependence': {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}': { + # 'model': 'polynomial', + # 'enthalpy-coefficients': {'class': 'np_array', 'object': [(1,'J/mol'), (2,'J/mol'), (3,'J/mol')]}, + # 'entropy-coefficients': {'class': 'np_array', 'object': [(1,'J/(mol*K)'),(2,'J/(mol*K)'),(3,'J/(mol*K)')]}} + # } } def test_make_wilhoit(self): From b9e8f129f26fbc98dc86afbddf73322379622ad1 Mon Sep 17 00:00:00 2001 From: Bjarne Kreitz Date: Tue, 14 Oct 2025 10:14:39 -0400 Subject: [PATCH 38/49] fix solverSurfaceTest.py by adding units --- test/rmgpy/solver/solverSurfaceTest.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/rmgpy/solver/solverSurfaceTest.py b/test/rmgpy/solver/solverSurfaceTest.py index 42fb9c81404..d4af80264f4 100644 --- a/test/rmgpy/solver/solverSurfaceTest.py +++ b/test/rmgpy/solver/solverSurfaceTest.py @@ -802,7 +802,7 @@ def test_solve_h2_thermo_coverage_dependence(self): Cpdata=([1.50, 2.58, 3.40, 4.00, 4.73, 5.13, 5.57], "cal/(mol*K)"), H298=(-11.26, "kcal/mol"), S298=(0.44, "cal/(mol*K)"), - thermo_coverage_dependence={"1 H u0 p0 {2,S} \n 2 X u0 p0 {1,S}":{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,5,3]},} + thermo_coverage_dependence={"1 H u0 p0 {2,S} \n 2 X u0 p0 {1,S}":{'model':'polynomial', 'enthalpy-coefficients':[(1,'J/mol'),(2,'J/mol'),(3,'J/mol')], "entropy-coefficients":[(1,'J/(mol*K)'),(5,'J/(mol*K)'),(3,'J/(mol*K)')]},} ), ) @@ -1008,7 +1008,7 @@ def test_solve_ch3_thermo_coverage_dependence(self): Tmax=(2000, "K"), E0=(-39.1285, "kJ/mol"), thermo_coverage_dependence={"1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S} \n 2 H u0 p0 c0 {1,S} \n 3 H u0 p0 c0 {1,S} \n 4 H u0 p0 c0 {1,S} \n 5 X u0 p0 c0 {1,S}": - {'model':'polynomial', 'enthalpy-coefficients':[1e5,2,3], "entropy-coefficients":[1,5,3]},}, + {'model':'polynomial', 'enthalpy-coefficients':[(1e5,'J/mol'),(2,'J/mol'),(3,'J/mol')], "entropy-coefficients":[(1,'J/(mol*K)'),(5,'J/(mol*K)'),(3,'J/(mol*K)')]},}, comment="""Thermo library: surfaceThermoNi111""", ), ) From a6c1de969bc14cec45640eb2b22c944aec1a6701 Mon Sep 17 00:00:00 2001 From: Bjarne Kreitz Date: Tue, 14 Oct 2025 14:28:45 -0400 Subject: [PATCH 39/49] calculate rxns_free_energy_change via np.matmul --- rmgpy/solver/surface.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index 2fad5164f16..a34a8645c2f 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -489,7 +489,7 @@ cdef class SurfaceReactor(ReactionSystem): for matrix in self.thermo_coeff_matrix: sp_free_energy_correction = np.diag(np.dot(matrix, thermo_dep_coverage)).sum() free_energy_coverage_corrections.append(sp_free_energy_correction) - rxns_free_energy_change = np.diag(np.dot(self.stoi_matrix, np.transpose(np.array([free_energy_coverage_corrections])))) + 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 From fa1196e4b91c5354d3dece9ec562a8049fde85d7 Mon Sep 17 00:00:00 2001 From: Bjarne Kreitz Date: Tue, 14 Oct 2025 14:29:49 -0400 Subject: [PATCH 40/49] remove Chaos bugfix that the coverage kinetics cannot be written to chemkin --- rmgpy/chemkin.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index 8b263fd0186..1988419057d 100644 --- a/rmgpy/chemkin.pyx +++ b/rmgpy/chemkin.pyx @@ -1867,7 +1867,7 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals for species, cov_params in kinetics.coverage_dependence.items(): label = get_species_identifier(species) string += f' COV / {label:<41} ' - string += f"{cov_params['a'].value_si:<9.3g} {cov_params['m'].value_si:<9.3g} {cov_params['E'].value_si/4184.:<9.3f} /\n" + string += f"{cov_params['a']:<9.3g} {cov_params['m']:<9.3g} {cov_params['E'].value_si/4184.:<9.3f} /\n" if isinstance(kinetics, (_kinetics.ThirdBody, _kinetics.Lindemann, _kinetics.Troe)): # Write collider efficiencies From 261230f582940d2cd1c65ad52b0b0ee0b66f9f7a Mon Sep 17 00:00:00 2001 From: Bjarne Kreitz Date: Mon, 6 Apr 2026 12:12:50 -0400 Subject: [PATCH 41/49] remove cantera dependence --- rmgpy/rmg/main.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index 98a027fff12..ddef94eb033 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -49,7 +49,6 @@ import yaml from cantera import ck2yaml from scipy.optimize import brute -import cantera as ct import rmgpy.util as util from rmgpy import settings @@ -1235,14 +1234,13 @@ def execute(self, initialize=True, **kwargs): if self.thermo_coverage_dependence: # add thermo coverage dependence to Cantera files chem_yaml_path = os.path.join(self.output_directory, "cantera", "chem.yaml") - gas = ct.Solution(chem_yaml_path, "gas") - surf = ct.Interface(chem_yaml_path, "surface1", [gas]) with open(chem_yaml_path, 'r') as f: content = yaml.load(f, Loader=yaml.FullLoader) content['phases'][1]['reference-state-coverage'] = 0.11 content['phases'][1]['thermo'] = 'coverage-dependent-surface' - + cantera_names = [content["species"][i]['name'] for i in range(len(content["species"]))] + for s in self.reaction_model.core.species: if s.contains_surface_site() and s.thermo.thermo_coverage_dependence: for dep_sp, parameters in s.thermo.thermo_coverage_dependence.items(): @@ -1253,9 +1251,9 @@ def execute(self, initialize=True, **kwargs): parameters['enthalpy-coefficients'] = [value.value_si for value in parameters['enthalpy-coefficients']] parameters['entropy-coefficients'] = [value.value_si for value in parameters['entropy-coefficients']] try: - content["species"][gas.n_species+surf.species_index(sp.to_chemkin())]['coverage-dependencies'][sp.to_chemkin()] = parameters + content["species"][cantera_names.index(s.to_chemkin())]['coverage-dependencies'][sp.to_chemkin()] = parameters except KeyError: - content["species"][gas.n_species+surf.species_index(sp.to_chemkin())]['coverage-dependencies'] = {sp.to_chemkin(): parameters} + content["species"][cantera_names.index(s.to_chemkin())]['coverage-dependencies'] = {sp.to_chemkin(): parameters} annotated_yaml_path = os.path.join(self.output_directory, "cantera", "chem_annotated.yaml") with open(annotated_yaml_path, 'r') as f: From 634aa035ce12ca5d17f98d82b144c2d6221aa3b1 Mon Sep 17 00:00:00 2001 From: Bjarne Kreitz Date: Mon, 6 Apr 2026 15:17:45 -0400 Subject: [PATCH 42/49] remove np_list function --- rmgpy/thermo/nasa.pyx | 5 ++--- rmgpy/thermo/thermodata.pyx | 5 ++--- rmgpy/thermo/wilhoit.pyx | 5 ++--- rmgpy/util.py | 11 ----------- 4 files changed, 6 insertions(+), 20 deletions(-) diff --git a/rmgpy/thermo/nasa.pyx b/rmgpy/thermo/nasa.pyx index 0b04438654a..8a1df5fe9c5 100644 --- a/rmgpy/thermo/nasa.pyx +++ b/rmgpy/thermo/nasa.pyx @@ -34,7 +34,6 @@ cimport numpy as np from libc.math cimport log cimport rmgpy.constants as constants -from rmgpy.util import np_list import rmgpy.quantity as quantity ################################################################################ @@ -315,8 +314,8 @@ cdef class NASA(HeatCapacityModel): for species, parameters in value.items(): # just the polynomial model for now processed_parameters = {'model': parameters['model'], - 'enthalpy-coefficients': np_list([quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']]), - 'entropy-coefficients': np_list([quantity.Entropy(p) for p in parameters['entropy-coefficients']]), + 'enthalpy-coefficients': [quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']], + 'entropy-coefficients': [quantity.Entropy(p) for p in parameters['entropy-coefficients']], } self._thermo_coverage_dependence[species] = processed_parameters diff --git a/rmgpy/thermo/thermodata.pyx b/rmgpy/thermo/thermodata.pyx index 6797fa99e4d..61698cc174f 100644 --- a/rmgpy/thermo/thermodata.pyx +++ b/rmgpy/thermo/thermodata.pyx @@ -35,7 +35,6 @@ cimport numpy as np from libc.math cimport log import rmgpy.quantity as quantity -from rmgpy.util import np_list ################################################################################ @@ -129,8 +128,8 @@ cdef class ThermoData(HeatCapacityModel): for species, parameters in value.items(): # just the polynomial model for now processed_parameters = {'model': parameters['model'], - 'enthalpy-coefficients': np_list([quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']]), - 'entropy-coefficients': np_list([quantity.Entropy(p) for p in parameters['entropy-coefficients']]), + 'enthalpy-coefficients': [quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']], + 'entropy-coefficients': [quantity.Entropy(p) for p in parameters['entropy-coefficients']], } self._thermo_coverage_dependence[species] = processed_parameters diff --git a/rmgpy/thermo/wilhoit.pyx b/rmgpy/thermo/wilhoit.pyx index bcbf16d4abb..4ba3980b542 100644 --- a/rmgpy/thermo/wilhoit.pyx +++ b/rmgpy/thermo/wilhoit.pyx @@ -35,7 +35,6 @@ from libc.math cimport sqrt, log cimport rmgpy.constants as constants import rmgpy.quantity as quantity -from rmgpy.util import np_list # Prior to numpy 1.14, `numpy.linalg.lstsq` does not accept None as a value RCOND = -1 if int(np.__version__.split('.')[1]) < 14 else None @@ -151,8 +150,8 @@ cdef class Wilhoit(HeatCapacityModel): for species, parameters in value.items(): # just the polynomial model for now processed_parameters = {'model': parameters['model'], - 'enthalpy-coefficients': np_list([quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']]), - 'entropy-coefficients': np_list([quantity.Entropy(p) for p in parameters['entropy-coefficients']]), + 'enthalpy-coefficients': [quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']], + 'entropy-coefficients': [quantity.Entropy(p) for p in parameters['entropy-coefficients']], } self._thermo_coverage_dependence[species] = processed_parameters diff --git a/rmgpy/util.py b/rmgpy/util.py index 365833379ea..15769e02e61 100644 --- a/rmgpy/util.py +++ b/rmgpy/util.py @@ -33,7 +33,6 @@ import shutil import time from functools import wraps -import numpy as np class Subject(object): @@ -239,13 +238,3 @@ def as_list(item, default=None): return default else: return [item] - -class np_list(np.ndarray): - """ - A subclass of numpy.ndarray which rendered as a list when printed. - """ - def __new__(cls, input_array): - obj = np.asarray(input_array).view(cls) - return obj - def __repr__(self): - return str(self.tolist()) From c15694c9046b6d08b78c82850ad7eec849f58be2 Mon Sep 17 00:00:00 2001 From: Richard West Date: Wed, 8 Apr 2026 11:02:31 -0400 Subject: [PATCH 43/49] Optimization: avoiding exponentials and stacks. --- rmgpy/solver/surface.pyx | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index a34a8645c2f..345138c72b2 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -483,8 +483,16 @@ cdef class SurfaceReactor(ReactionSystem): else: surface_site_fraction = 0 coverages.append(surface_site_fraction) - coverages = np.array(coverages) - thermo_dep_coverage = np.stack([coverages, coverages**2, coverages**3, -self.T.value_si*coverages, -self.T.value_si*coverages**2, -self.T.value_si*coverages**3]) + coverages = np.array(coverages, dtype=np.float64) + 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 = [] for matrix in self.thermo_coeff_matrix: sp_free_energy_correction = np.diag(np.dot(matrix, thermo_dep_coverage)).sum() From e185c4e7a7b7ca3ba563c8fd4d96bae43d071572 Mon Sep 17 00:00:00 2001 From: Richard West Date: Wed, 8 Apr 2026 11:05:38 -0400 Subject: [PATCH 44/49] optimization: removing loops, and making cython declarations --- rmgpy/solver/surface.pyx | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index 345138c72b2..e19cd3bbb4a 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -433,6 +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 @@ -476,14 +480,7 @@ cdef class SurfaceReactor(ReactionSystem): # Thermodynamic coverage dependence if self.thermo_coverage_dependence: - coverages = [] - for i in range(len(N)): - if species_on_surface[i]: - surface_site_fraction = N[i] / total_sites - else: - surface_site_fraction = 0 - coverages.append(surface_site_fraction) - coverages = np.array(coverages, dtype=np.float64) + 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) @@ -493,10 +490,9 @@ cdef class SurfaceReactor(ReactionSystem): 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 = [] - for matrix in self.thermo_coeff_matrix: - sp_free_energy_correction = np.diag(np.dot(matrix, thermo_dep_coverage)).sum() - free_energy_coverage_corrections.append(sp_free_energy_correction) + 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)) From 1386801e50dd98c11b25c76a5f63efd6958f6d7a Mon Sep 17 00:00:00 2001 From: Richard West Date: Wed, 8 Apr 2026 11:07:37 -0400 Subject: [PATCH 45/49] Don't save empty dictionary when thermo coverage is None Previously the as_dict and the __repr__ would both return thermo_coverage_dependence={} for absolutely everything. Now if it starts as None, it stays None. --- rmgpy/thermo/nasa.pyx | 18 ++++++++++-------- rmgpy/thermo/thermodata.pyx | 16 +++++++++------- rmgpy/thermo/wilhoit.pyx | 16 +++++++++------- 3 files changed, 28 insertions(+), 22 deletions(-) diff --git a/rmgpy/thermo/nasa.pyx b/rmgpy/thermo/nasa.pyx index 8a1df5fe9c5..0ee505e4167 100644 --- a/rmgpy/thermo/nasa.pyx +++ b/rmgpy/thermo/nasa.pyx @@ -309,16 +309,18 @@ cdef class NASA(HeatCapacityModel): def __get__(self): return self._thermo_coverage_dependence def __set__(self, value): - self._thermo_coverage_dependence = {} if value: - for species, parameters in value.items(): + self._thermo_coverage_dependence = {} + for species, parameters in value.items(): # just the polynomial model for now - processed_parameters = {'model': parameters['model'], - 'enthalpy-coefficients': [quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']], - 'entropy-coefficients': [quantity.Entropy(p) for p in parameters['entropy-coefficients']], - } - self._thermo_coverage_dependence[species] = processed_parameters - + processed_parameters = {'model': parameters['model'], + 'enthalpy-coefficients': [quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']], + 'entropy-coefficients': [quantity.Entropy(p) for p in parameters['entropy-coefficients']], + } + self._thermo_coverage_dependence[species] = processed_parameters + else: + self._thermo_coverage_dependence = None + cpdef double get_heat_capacity(self, double T) except -1000000000: """ Return the constant-pressure heat capacity diff --git a/rmgpy/thermo/thermodata.pyx b/rmgpy/thermo/thermodata.pyx index 61698cc174f..0b2afb3a69d 100644 --- a/rmgpy/thermo/thermodata.pyx +++ b/rmgpy/thermo/thermodata.pyx @@ -123,15 +123,17 @@ cdef class ThermoData(HeatCapacityModel): def __get__(self): return self._thermo_coverage_dependence def __set__(self, value): - self._thermo_coverage_dependence = {} if value: - for species, parameters in value.items(): + self._thermo_coverage_dependence = {} + for species, parameters in value.items(): # just the polynomial model for now - processed_parameters = {'model': parameters['model'], - 'enthalpy-coefficients': [quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']], - 'entropy-coefficients': [quantity.Entropy(p) for p in parameters['entropy-coefficients']], - } - self._thermo_coverage_dependence[species] = processed_parameters + processed_parameters = {'model': parameters['model'], + 'enthalpy-coefficients': [quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']], + 'entropy-coefficients': [quantity.Entropy(p) for p in parameters['entropy-coefficients']], + } + self._thermo_coverage_dependence[species] = processed_parameters + else: + self._thermo_coverage_dependence = None @cython.boundscheck(False) @cython.wraparound(False) diff --git a/rmgpy/thermo/wilhoit.pyx b/rmgpy/thermo/wilhoit.pyx index 4ba3980b542..e4666a8ecfb 100644 --- a/rmgpy/thermo/wilhoit.pyx +++ b/rmgpy/thermo/wilhoit.pyx @@ -145,15 +145,17 @@ cdef class Wilhoit(HeatCapacityModel): def __get__(self): return self._thermo_coverage_dependence def __set__(self, value): - self._thermo_coverage_dependence = {} if value: - for species, parameters in value.items(): + self._thermo_coverage_dependence = {} + for species, parameters in value.items(): # just the polynomial model for now - processed_parameters = {'model': parameters['model'], - 'enthalpy-coefficients': [quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']], - 'entropy-coefficients': [quantity.Entropy(p) for p in parameters['entropy-coefficients']], - } - self._thermo_coverage_dependence[species] = processed_parameters + processed_parameters = {'model': parameters['model'], + 'enthalpy-coefficients': [quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']], + 'entropy-coefficients': [quantity.Entropy(p) for p in parameters['entropy-coefficients']], + } + self._thermo_coverage_dependence[species] = processed_parameters + else: + self._thermo_coverage_dependence = None cpdef double get_heat_capacity(self, double T) except -1000000000: """ From 84b17499bb33b0eae7ea5ccde3629ea92b7d9fae Mon Sep 17 00:00:00 2001 From: Richard West Date: Wed, 8 Apr 2026 11:08:27 -0400 Subject: [PATCH 46/49] Small optimization and simplification. --- rmgpy/rmg/main.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index ddef94eb033..e38218a6ce5 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -1239,7 +1239,7 @@ def execute(self, initialize=True, **kwargs): content['phases'][1]['reference-state-coverage'] = 0.11 content['phases'][1]['thermo'] = 'coverage-dependent-surface' - cantera_names = [content["species"][i]['name'] for i in range(len(content["species"]))] + cantera_names = [s['name'] for s in content["species"]] for s in self.reaction_model.core.species: if s.contains_surface_site() and s.thermo.thermo_coverage_dependence: From ea374fa9acf3e8bac17e52870e9ba7f675f87d9d Mon Sep 17 00:00:00 2001 From: Richard West Date: Wed, 8 Apr 2026 12:22:29 -0400 Subject: [PATCH 47/49] update tests for coverage-dependence data as_dict The on-disk/YAML format for enthalpy-coefficients/entropy-coefficients changed: Old: {'class': 'np_array', 'object': [, ...]} New: [{'class': 'ScalarQuantity', 'value': ..., 'units': '...'}, ...] --- test/rmgpy/thermo/nasaTest.py | 9 +++++---- test/rmgpy/thermo/wilhoitTest.py | 8 ++++---- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/test/rmgpy/thermo/nasaTest.py b/test/rmgpy/thermo/nasaTest.py index 138ea59d3f7..e7e59de601b 100644 --- a/test/rmgpy/thermo/nasaTest.py +++ b/test/rmgpy/thermo/nasaTest.py @@ -378,10 +378,10 @@ def test_nasa_as_dict_full(self): assert nasa_dict["thermo_coverage_dependence"].keys() == self.thermo_coverage_dependence.keys() sp_name = list(self.thermo_coverage_dependence.keys())[0] assert nasa_dict['thermo_coverage_dependence'][sp_name]['model'] == self.thermo_coverage_dependence[sp_name]['model'] - enthalpy_list = nasa_dict['thermo_coverage_dependence'][sp_name]['enthalpy-coefficients']['object'] - assert [(int(coeff.value), str(coeff.units)) for coeff in enthalpy_list] == self.thermo_coverage_dependence[sp_name]['enthalpy-coefficients'] - entropy_list = nasa_dict['thermo_coverage_dependence'][sp_name]['entropy-coefficients']['object'] - assert [(int(coeff.value), str(coeff.units)) for coeff in entropy_list] == self.thermo_coverage_dependence[sp_name]['entropy-coefficients'] + enthalpy_list = nasa_dict['thermo_coverage_dependence'][sp_name]['enthalpy-coefficients'] + assert [(int(coeff['value']), str(coeff['units'])) for coeff in enthalpy_list] == self.thermo_coverage_dependence[sp_name]['enthalpy-coefficients'] + entropy_list = nasa_dict['thermo_coverage_dependence'][sp_name]['entropy-coefficients'] + assert [(int(coeff['value']), str(coeff['units'])) for coeff in entropy_list] == self.thermo_coverage_dependence[sp_name]['entropy-coefficients'] assert nasa_dict["comment"] == self.comment assert tuple(nasa_dict["polynomials"]["polynomial1"]["coeffs"]["object"]) == tuple(self.coeffs_low) assert tuple(nasa_dict["polynomials"]["polynomial2"]["coeffs"]["object"]) == tuple(self.coeffs_high) @@ -403,6 +403,7 @@ def test_nasa_as_dict_minimal(self): assert "CpInf" not in keys assert "label" not in keys assert "comment" not in keys + assert "thermo_coverage_dependence" not in keys def test_nasa_polynomial_as_dict(self): """ diff --git a/test/rmgpy/thermo/wilhoitTest.py b/test/rmgpy/thermo/wilhoitTest.py index 92fe4737480..3b8844c184e 100644 --- a/test/rmgpy/thermo/wilhoitTest.py +++ b/test/rmgpy/thermo/wilhoitTest.py @@ -432,10 +432,10 @@ def test_wilhoit_as_dict(self): assert wilhoit_dict["thermo_coverage_dependence"].keys() == self.thermo_coverage_dependence.keys() sp_name = list(self.thermo_coverage_dependence.keys())[0] assert wilhoit_dict['thermo_coverage_dependence'][sp_name]['model'] == self.thermo_coverage_dependence[sp_name]['model'] - enthalpy_list = wilhoit_dict['thermo_coverage_dependence'][sp_name]['enthalpy-coefficients']['object'] - assert [(int(coeff.value), str(coeff.units)) for coeff in enthalpy_list] == self.thermo_coverage_dependence[sp_name]['enthalpy-coefficients'] - entropy_list = wilhoit_dict['thermo_coverage_dependence'][sp_name]['entropy-coefficients']['object'] - assert [(int(coeff.value), str(coeff.units)) for coeff in entropy_list] == self.thermo_coverage_dependence[sp_name]['entropy-coefficients'] + enthalpy_list = wilhoit_dict['thermo_coverage_dependence'][sp_name]['enthalpy-coefficients'] + assert [(int(coeff['value']), str(coeff['units'])) for coeff in enthalpy_list] == self.thermo_coverage_dependence[sp_name]['enthalpy-coefficients'] + entropy_list = wilhoit_dict['thermo_coverage_dependence'][sp_name]['entropy-coefficients'] + assert [(int(coeff['value']), str(coeff['units'])) for coeff in entropy_list] == self.thermo_coverage_dependence[sp_name]['entropy-coefficients'] wilhoit_dict = {k: wilhoit_dict[k] for k in wilhoit_dict.keys() - {'thermo_coverage_dependence'}} assert wilhoit_dict == { From 589fe738f63ad4a003bd31f65f2c9e9c52b71c7e Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 9 Apr 2026 11:38:03 -0400 Subject: [PATCH 48/49] Fix YAML dumping style for consistency in output files. The default_flow_style=None will let it use a single line for things like [1.3, 4.0, 4.45]. Which should look more like we expected. It's a bit clunky that we have to read the yaml, modify the dicts, and dump the yaml, but this should make the changes less obnoxious. --- rmgpy/rmg/main.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index e38218a6ce5..dc7100a4367 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -1263,10 +1263,10 @@ def execute(self, initialize=True, **kwargs): annotated_content['species'] = content['species'] with open(chem_yaml_path, 'w') as output_f: - yaml.dump(content, output_f, sort_keys=False) + yaml.dump(content, output_f, sort_keys=False, default_flow_style=None) with open(annotated_yaml_path, 'w') as output_f: - yaml.dump(annotated_content, output_f, sort_keys=False) + yaml.dump(annotated_content, output_f, sort_keys=False, default_flow_style=None) else: # gas phase only self.generate_cantera_files(os.path.join(self.output_directory, "chemkin", "chem.inp")) From c9e20b159173d60aaefd3362bef0db72e0dce893 Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 9 Apr 2026 15:41:53 -0400 Subject: [PATCH 49/49] Changing coverage-dependent thermo cantera yaml output. Now it just modifies the part of the yaml file that has changed, instead of reading and writing the whole thing. Will probably still do some refactoring, but making an intermediate commit to check it's working so far. Second commit: (squashed) A refactor of the coverage-dependence-thermo to cantera yaml code. Moved the _add_coverage_dependence_to_cantera_yaml function, and made it pass strings instead of dictionaries. --- rmgpy/rmg/main.py | 96 +++++++++++++++++++++++++++++++---------------- 1 file changed, 64 insertions(+), 32 deletions(-) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index dc7100a4367..3ce2ce41202 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -37,6 +37,7 @@ import logging import marshal import os +import re import resource import shutil import sys @@ -1223,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")), @@ -1231,43 +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: - # add thermo coverage dependence to Cantera files - chem_yaml_path = os.path.join(self.output_directory, "cantera", "chem.yaml") - with open(chem_yaml_path, 'r') as f: - content = yaml.load(f, Loader=yaml.FullLoader) - - content['phases'][1]['reference-state-coverage'] = 0.11 - content['phases'][1]['thermo'] = 'coverage-dependent-surface' - cantera_names = [s['name'] for s in content["species"]] + 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: - for dep_sp, parameters in s.thermo.thermo_coverage_dependence.items(): - mol = Molecule().from_adjacency_list(dep_sp) + 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): - parameters['units'] = {'energy':'J', 'quantity':'mol'} - parameters['enthalpy-coefficients'] = [value.value_si for value in parameters['enthalpy-coefficients']] - parameters['entropy-coefficients'] = [value.value_si for value in parameters['entropy-coefficients']] - try: - content["species"][cantera_names.index(s.to_chemkin())]['coverage-dependencies'][sp.to_chemkin()] = parameters - except KeyError: - content["species"][cantera_names.index(s.to_chemkin())]['coverage-dependencies'] = {sp.to_chemkin(): parameters} - - annotated_yaml_path = os.path.join(self.output_directory, "cantera", "chem_annotated.yaml") - with open(annotated_yaml_path, 'r') as f: - annotated_content = yaml.load(f, Loader=yaml.FullLoader) - - annotated_content['phases'] = content['phases'] - annotated_content['species'] = content['species'] - - with open(chem_yaml_path, 'w') as output_f: - yaml.dump(content, output_f, sort_keys=False, default_flow_style=None) - - with open(annotated_yaml_path, 'w') as output_f: - yaml.dump(annotated_content, output_f, sort_keys=False, default_flow_style=None) - + 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")) @@ -2431,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): """