From 166c13acdc8894e1ab63bd8acd197ba81103f4a5 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Sun, 10 Mar 2024 16:49:42 -0400 Subject: [PATCH 01/30] 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 105a1d68f7..3c20de6ab9 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 e7f6b082b5ed4daa9169830756355894372a7409 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 11 Mar 2024 13:26:58 -0400 Subject: [PATCH 02/30] 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 e04cd82e2e..9d0341d2a6 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -105,7 +105,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. @@ -152,6 +153,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): """ @@ -1062,7 +1064,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 27ab4d56b30decf01adf4050d867cda0c5fc2f69 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 11 Mar 2024 13:27:21 -0400 Subject: [PATCH 03/30] 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 3c20de6ab9..eb9f2b0996 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 95904258f3cca021dc33f7cb0f50a2f62fe3c30f Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 11 Mar 2024 15:30:01 -0400 Subject: [PATCH 04/30] 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 eb9f2b0996..54151c4993 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 dcb95ca45051bf6fa68b5166b38df63389f5610a Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 11 Mar 2024 20:43:28 -0400 Subject: [PATCH 05/30] 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 e484f3036f..832183ac75 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 5248e79a56..99c7237789 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 20a1061b18..2dc974b016 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 6cd68ea3b1deadbdb86f36307f2fea57ae8d2cff Mon Sep 17 00:00:00 2001 From: 12Chao Date: Tue, 12 Mar 2024 14:47:25 -0400 Subject: [PATCH 06/30] 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 b462301f72..fb221419da 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 d9052db0b4..0a1d08b12f 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 576c8a2e67..e2cd2419e4 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 54bbaa3d4eb7cfde18145d7a3f6d8e5e26af006f Mon Sep 17 00:00:00 2001 From: 12Chao Date: Tue, 12 Mar 2024 14:48:37 -0400 Subject: [PATCH 07/30] 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 832183ac75..8b2a65ae94 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 99c7237789..e38ae87b51 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 2dc974b016..1529966636 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 adf805241e711e5a181722834c50e4e99750f8a6 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Tue, 12 Mar 2024 16:12:12 -0400 Subject: [PATCH 08/30] 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 fb221419da..43f629e60c 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 8b2a65ae94..d7af53acd0 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 542857bda820ba92db28b8ad8f3d8140ca1789c6 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Wed, 13 Mar 2024 14:53:34 -0400 Subject: [PATCH 09/30] 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 43f629e60c..6922c51e54 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 d7af53acd0..4014cdfc1b 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 3de057a193d0307c99922575b6d8ab6b0a6d04f7 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Wed, 13 Mar 2024 16:01:00 -0400 Subject: [PATCH 10/30] 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 94637ce8f6..d35f3de4c7 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 f08b28e0da..8945bef5f4 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 e3a633c2da..cca336e133 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 4b74fd72a7..15d76d34fa 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 7d4ede92a192c412766e7d7ef7c2c0b28e069054 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Wed, 13 Mar 2024 16:11:59 -0400 Subject: [PATCH 11/30] 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 d35f3de4c7..ef2f7265e3 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 02be9566c0109367582c1c870a79876f667fc00d Mon Sep 17 00:00:00 2001 From: 12Chao Date: Wed, 13 Mar 2024 19:19:42 -0400 Subject: [PATCH 12/30] 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 54151c4993..4b15a0240f 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 b15d065566353a47aad70a8f790ce78478f08ca5 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Fri, 15 Mar 2024 17:58:21 -0400 Subject: [PATCH 13/30] 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 4b15a0240f..ab334546b8 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 cfa2babb431e7198ba9d333139a12bbc76315780 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Fri, 15 Mar 2024 21:25:26 -0400 Subject: [PATCH 14/30] 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 ab334546b8..b943df8263 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 5a2c8d80f20023337cc810245e414353e1545758 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 18 Mar 2024 19:24:25 -0400 Subject: [PATCH 15/30] 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 b943df8263..cddd981d27 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 7f2be04054229d57373940c9caeee73399a65628 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Wed, 20 Mar 2024 12:05:02 -0400 Subject: [PATCH 16/30] 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 36c08b4ad8..dc9ee25b3f 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 4fa77159629ecceb1acdc3f468a491bd84ebe55f Mon Sep 17 00:00:00 2001 From: 12Chao Date: Sun, 24 Mar 2024 19:59:55 -0400 Subject: [PATCH 17/30] 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 cddd981d27..20b01311bf 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 da3f9e1895b5bddba5fa45990c9c88665968e101 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 25 Mar 2024 12:59:18 -0400 Subject: [PATCH 18/30] 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 20b01311bf..becdd844f4 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 9f26e2c66b80aa3ead85ce20aad9b80d3326e4f6 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 25 Mar 2024 13:00:30 -0400 Subject: [PATCH 19/30] 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 dc9ee25b3f..42fb9c8140 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 6ca15a190dee90a795a0095aad9e477a07446a2c Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 25 Mar 2024 14:12:03 -0400 Subject: [PATCH 20/30] 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 4014cdfc1b..d749f986cb 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 e38ae87b51..9274e088f9 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 1529966636..8801511491 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 1ab80e37201249c4e4c01a8206b1fdbe42b69502 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 25 Mar 2024 17:47:05 -0400 Subject: [PATCH 21/30] 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 15769e02e6..365833379e 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 80974dcfcdd871d5a6a851bb7a7780f3df737527 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 25 Mar 2024 17:52:59 -0400 Subject: [PATCH 22/30] 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 d749f986cb..e270400ae5 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 9274e088f9..aed7d18f89 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 8801511491..163b39e187 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 fcc050be832b47bffa23561401f7bc9e05e0c6d1 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 25 Mar 2024 17:53:56 -0400 Subject: [PATCH 23/30] 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 8945bef5f4..f96df94cfe 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 cca336e133..3c6b39e7e0 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 15d76d34fa..6a350d7833 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 3f190d3ec03a7f77e4e89d4e2eacba7a1dc330f6 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Mon, 25 Mar 2024 18:01:32 -0400 Subject: [PATCH 24/30] 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 becdd844f4..1c60273406 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 d5f5349ae763b93c519a7835195f84244ca783e3 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Thu, 28 Mar 2024 14:16:08 -0400 Subject: [PATCH 25/30] 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 8df771980d..805e97230a 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -50,6 +50,7 @@ import yaml from cantera import ck2yaml from scipy.optimize import brute +import cantera as ct import rmgpy.util as util from rmgpy.rmg.model import Species, CoreEdgeReactionModel @@ -1179,6 +1180,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 1900513448b9dce9948142e5afec2a74c4dbc530 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Fri, 29 Mar 2024 14:29:47 -0400 Subject: [PATCH 26/30] 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 805e97230a..bb9fa44e78 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -1199,9 +1199,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 984bc6cd2b7c1c8aabd4a5b61c81da626e5dd4ce Mon Sep 17 00:00:00 2001 From: 12Chao Date: Fri, 29 Mar 2024 21:06:43 -0400 Subject: [PATCH 27/30] 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 bb9fa44e78..a163ae16a1 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -1197,11 +1197,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 21ad0b03996fb257ca1dd94a37becde2922f8225 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Tue, 2 Apr 2024 12:50:51 -0400 Subject: [PATCH 28/30] Add thermo_coverag_dependence as an instance for class RMG --- rmgpy/rmg/main.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index a163ae16a1..630fcdd6c4 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -191,6 +191,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 @@ -275,6 +276,7 @@ def load_input(self, path=None): 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 self.reaction_model.verbose_comments = self.verbose_comments self.reaction_model.save_edge_species = self.save_edge_species From e41093ec1b915a59e12978f8c8f731cae2a64afa Mon Sep 17 00:00:00 2001 From: 12Chao Date: Wed, 17 Jul 2024 20:40:46 -0400 Subject: [PATCH 29/30] 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 d28ee47ce0..f225f078d3 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 426af6a0b1dc6996952ce63d76f557bd48459784 Mon Sep 17 00:00:00 2001 From: 12Chao Date: Sat, 19 Oct 2024 10:18:17 -0400 Subject: [PATCH 30/30] 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 630fcdd6c4..1e575bf800 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -1200,8 +1200,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 1c60273406..2fad5164f1 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 e270400ae5..0b04438654 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 aed7d18f89..6797fa99e4 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 163b39e187..3342f64c48 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 f96df94cfe..23255082f9 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)