From 67f934da2aebf8fa07d1d1a5049270fdff9e1316 Mon Sep 17 00:00:00 2001 From: Sevy Harris Date: Thu, 19 Oct 2023 14:44:59 -0400 Subject: [PATCH 1/6] Add surface mechanisms to regression tool --- rmgpy/tools/canteramodel.py | 118 ++++++++++++++++++++++----- rmgpy/tools/observablesregression.py | 116 +++++++++++++++++++++----- rmgpy/tools/regression.py | 19 ++++- 3 files changed, 211 insertions(+), 42 deletions(-) diff --git a/rmgpy/tools/canteramodel.py b/rmgpy/tools/canteramodel.py index 99aed79de1..15d26c2947 100644 --- a/rmgpy/tools/canteramodel.py +++ b/rmgpy/tools/canteramodel.py @@ -28,6 +28,7 @@ ############################################################################### import os.path +import logging import cantera as ct import numpy as np @@ -63,18 +64,44 @@ class CanteraCondition(object): """ - def __init__(self, reactor_type, reaction_time, mol_frac, T0=None, P0=None, V0=None): + def __init__(self, reactor_type, reaction_time, mol_frac, surface_mol_frac=None, T0=None, P0=None, V0=None): self.reactor_type = reactor_type self.reaction_time = Quantity(reaction_time) # Normalize initialMolFrac if not already done: if sum(mol_frac.values()) != 1.00: + logging.warning('Initial mole fractions do not sum to one; normalizing.') + logging.info('') + logging.info('Original composition:') + for spec, molfrac in mol_frac.items(): + logging.info(f'{spec} = {molfrac}') total = sum(mol_frac.values()) + logging.info('') + logging.info('Normalized mole fractions:') for species, value in mol_frac.items(): mol_frac[species] = value / total + logging.info(f'{species} = {mol_frac[species]}') + logging.info('') self.mol_frac = mol_frac + if surface_mol_frac: + # normalize surface mol fractions + if sum(surface_mol_frac.values()) != 1.00: + logging.warning('Initial surface mole fractions do not sum to one; normalizing.') + logging.info('') + logging.info('Original composition:') + for spec, molfrac in surface_mol_frac.items(): + logging.info(f'{spec} = {molfrac}') + logging.info('') + logging.info('Normalized mole fractions:') + total = sum(surface_mol_frac.values()) + for species, value in surface_mol_frac.items(): + surface_mol_frac[species] = value / total + logging.info(f'{species} = {surface_mol_frac[species]}') + logging.info('') + self.surface_mol_frac = surface_mol_frac + # Check to see that one of the three attributes T0, P0, and V0 is less unspecified props = [T0, P0, V0] total = 0 @@ -121,7 +148,7 @@ def __str__(self): return string -def generate_cantera_conditions(reactor_type_list, reaction_time_list, mol_frac_list, Tlist=None, Plist=None, Vlist=None): +def generate_cantera_conditions(reactor_type_list, reaction_time_list, mol_frac_list, surface_mol_frac_list=None, Tlist=None, Plist=None, Vlist=None): """ Creates a list of cantera conditions from from the arguments provided. @@ -137,6 +164,8 @@ def generate_cantera_conditions(reactor_type_list, reaction_time_list, mol_frac_ `reaction_time_list` A tuple object giving the ([list of reaction times], units) `mol_frac_list` A list of molfrac dictionaries with species object keys and mole fraction values + `surface_mol_frac_list` A list of molfrac dictionaries with surface species object keys + and mole fraction values To specify the system for an ideal gas, you must define 2 of the following 3 parameters: `T0List` A tuple giving the ([list of initial temperatures], units) 'P0List' A tuple giving the ([list of initial pressures], units) @@ -162,30 +191,35 @@ def generate_cantera_conditions(reactor_type_list, reaction_time_list, mol_frac_ for i in range(len(reaction_time_list.value))] conditions = [] + if surface_mol_frac_list is None: + surface_mol_frac_list = [] # initialize here to avoid mutable default argument if Tlist is None: for reactor_type in reactor_type_list: for reaction_time in reaction_time_list: for mol_frac in mol_frac_list: - for P in Plist: - for V in Vlist: - conditions.append(CanteraCondition(reactor_type, reaction_time, mol_frac, P0=P, V0=V)) + for surface_mol_frac in surface_mol_frac_list: + for P in Plist: + for V in Vlist: + conditions.append(CanteraCondition(reactor_type, reaction_time, mol_frac, surface_mol_frac=surface_mol_frac, P0=P, V0=V)) elif Plist is None: for reactor_type in reactor_type_list: for reaction_time in reaction_time_list: for mol_frac in mol_frac_list: - for T in Tlist: - for V in Vlist: - conditions.append(CanteraCondition(reactor_type, reaction_time, mol_frac, T0=T, V0=V)) + for surface_mol_frac in surface_mol_frac_list: + for T in Tlist: + for V in Vlist: + conditions.append(CanteraCondition(reactor_type, reaction_time, mol_frac, surface_mol_frac=surface_mol_frac, T0=T, V0=V)) elif Vlist is None: for reactor_type in reactor_type_list: for reaction_time in reaction_time_list: for mol_frac in mol_frac_list: - for T in Tlist: - for P in Plist: - conditions.append(CanteraCondition(reactor_type, reaction_time, mol_frac, T0=T, P0=P)) + for surface_mol_frac in surface_mol_frac_list: + for T in Tlist: + for P in Plist: + conditions.append(CanteraCondition(reactor_type, reaction_time, mol_frac, surface_mol_frac=surface_mol_frac, T0=T, P0=P)) else: raise Exception("Cantera conditions must leave one of T0, P0, and V0 state variables unspecified") @@ -198,7 +232,7 @@ class Cantera(object): """ def __init__(self, species_list=None, reaction_list=None, canteraFile='', output_directory='', conditions=None, - sensitive_species=None, thermo_SA=False): + sensitive_species=None, thermo_SA=False, surface_species_list=None, surface_reaction_list=None): """ `species_list`: list of RMG species objects `reaction_list`: list of RMG reaction objects @@ -208,16 +242,28 @@ def __init__(self, species_list=None, reaction_list=None, canteraFile='', output `sensitive_species`: a list of RMG species objects for conductng sensitivity analysis on `thermo_SA`: a boolean indicating whether or not to run thermo SA. By default, if sensitive_species is given, only kinetic_SA will be calculated and it must be additionally specified to perform thermo SA. + `surface_species_list`: list of RMG species objects contianing the surface species. This is necessary for all surface mechanisms to initialize the starting mole fractions. + `surface_reaction_list`: list of RMG reaction objects containing the surface reactions. + For surface mechanisms, either canteraFile must be provided, or load_chemkin_model() must be called later on to create the Cantera file with the surface mechanism. """ self.species_list = species_list self.reaction_list = reaction_list self.reaction_map = {} self.model = ct.Solution(canteraFile) if canteraFile else None + self.surface = bool(surface_species_list or surface_reaction_list) + self.surface_species_list = surface_species_list + self.surface_reaction_list = surface_reaction_list self.output_directory = output_directory if output_directory else os.getcwd() self.conditions = conditions if conditions else [] self.sensitive_species = sensitive_species if sensitive_species else [] self.thermo_SA = thermo_SA + if self.surface: + assert surface_species_list, "Surface species list must be specified to run a surface simulation" + if canteraFile: + self.surface = ct.Interface(canteraFile, 'surface1') + self.model = self.surface.adjacent['gas'] + # Make output directory if it does not yet exist: if not os.path.exists(self.output_directory): try: @@ -225,7 +271,7 @@ def __init__(self, species_list=None, reaction_list=None, canteraFile='', output except: raise Exception('Cantera output directory could not be created.') - def generate_conditions(self, reactor_type_list, reaction_time_list, mol_frac_list, Tlist=None, Plist=None, Vlist=None): + def generate_conditions(self, reactor_type_list, reaction_time_list, mol_frac_list, surface_mol_frac_list=None, Tlist=None, Plist=None, Vlist=None): """ This saves all the reaction conditions into the Cantera class. ======================= ==================================================== @@ -240,19 +286,19 @@ def generate_conditions(self, reactor_type_list, reaction_time_list, mol_frac_li `reaction_time_list` A tuple object giving the ([list of reaction times], units) `mol_frac_list` A list of molfrac dictionaries with species object keys and mole fraction values + `surface_mol_frac_list` A list of molfrac dictionaries with surface species object keys and mole fraction values To specify the system for an ideal gas, you must define 2 of the following 3 parameters: `T0List` A tuple giving the ([list of initial temperatures], units) 'P0List' A tuple giving the ([list of initial pressures], units) 'V0List' A tuple giving the ([list of initial specific volumes], units) """ - self.conditions = generate_cantera_conditions(reactor_type_list, reaction_time_list, mol_frac_list, Tlist, Plist, Vlist) + self.conditions = generate_cantera_conditions(reactor_type_list, reaction_time_list, mol_frac_list, surface_mol_frac_list, Tlist, Plist, Vlist) def load_model(self): """ Load a cantera Solution model from the job's own species_list and reaction_list attributes """ - ct_species = [spec.to_cantera(use_chemkin_identifier=True) for spec in self.species_list] self.reaction_map = {} @@ -292,6 +338,7 @@ def load_chemkin_model(self, chemkin_file, transport_file=None, **kwargs): Convert a chemkin mechanism chem.inp file to a cantera mechanism file chem.yaml and save it in the output_directory Then load it into self.model + Note that if this is a surface mechanism, the calling function much include surface_file as a keyword argument for the parser """ from cantera import ck2yaml @@ -301,9 +348,14 @@ def load_chemkin_model(self, chemkin_file, transport_file=None, **kwargs): if os.path.exists(out_name): os.remove(out_name) parser = ck2yaml.Parser() + parser.convert_mech(chemkin_file, transport_file=transport_file, out_name=out_name, **kwargs) self.model = ct.Solution(out_name) + if self.surface: + self.surface = ct.Interface(out_name, 'surface1') + self.model = self.surface.adjacent['gas'] + def modify_reaction_kinetics(self, rmg_reaction_index, rmg_reaction): """ Modify the corresponding cantera reaction's kinetics to match @@ -385,14 +437,23 @@ def simulate(self): species_names_list = [get_species_identifier(species) for species in self.species_list] inert_index_list = [self.species_list.index(species) for species in self.species_list if species.index == -1] + if self.surface: + surface_species_names_list = [get_species_identifier(species) for species in self.surface_species_list] + all_data = [] for condition in self.conditions: # First translate the mol_frac from species objects to species names new_mol_frac = {} - for key, value in condition.mol_frac.items(): - newkey = get_species_identifier(key) - new_mol_frac[newkey] = value + for rmg_species, mol_frac in condition.mol_frac.items(): + species_name = get_species_identifier(rmg_species) + new_mol_frac[species_name] = mol_frac + + if self.surface: + new_surface_mol_frac = {} + for rmg_species, mol_frac in condition.surface_mol_frac.items(): + species_name = get_species_identifier(rmg_species) + new_surface_mol_frac[species_name] = mol_frac # Set Cantera simulation conditions if condition.V0 is None: @@ -413,6 +474,11 @@ def simulate(self): else: raise Exception('Other types of reactor conditions are currently not supported') + # add the surface/gas adjacent thing... + if self.surface: + rsurf = ct.ReactorSurface(self.surface, cantera_reactor) + self.surface.TPX = condition.T0.value_si, condition.P0.value_si, new_surface_mol_frac + # Run this individual condition as a simulation cantera_simulation = ct.ReactorNet([cantera_reactor]) @@ -451,7 +517,12 @@ def simulate(self): times.append(cantera_simulation.time) temperature.append(cantera_reactor.T) pressure.append(cantera_reactor.thermo.P) - species_data.append(cantera_reactor.thermo[species_names_list].X) + + if self.surface: + species_data.append(np.concatenate((cantera_reactor.thermo[species_names_list].X, rsurf.kinetics.coverages))) + N_gas = len(cantera_reactor.thermo[species_names_list].X) + else: + species_data.append(cantera_reactor.thermo[species_names_list].X) if self.sensitive_species: # Cantera returns mass-based sensitivities rather than molar concentration or mole fraction based sensitivities. @@ -531,6 +602,15 @@ def simulate(self): index=species.index ) condition_data.append(species_generic_data) + if self.surface: + for index, species in enumerate(self.surface_species_list): + # Create generic data object that saves the surface species object into the species object. + species_generic_data = GenericData(label=surface_species_names_list[index], + species=species, + data=species_data[:, index + N_gas], + index=species.index + N_gas + ) + condition_data.append(species_generic_data) # save kinetic data as generic data object reaction_sensitivity_data = [] diff --git a/rmgpy/tools/observablesregression.py b/rmgpy/tools/observablesregression.py index 1cce9e1750..17819bb14e 100644 --- a/rmgpy/tools/observablesregression.py +++ b/rmgpy/tools/observablesregression.py @@ -112,32 +112,74 @@ def __init__(self, title='', old_dir='', new_dir='', observables=None, expt_data if os.path.exists(os.path.join(new_dir, 'tran.dat')): new_transport_path = os.path.join(new_dir, 'tran.dat') - # load the species and reactions from each model - old_species_list, old_reaction_list = load_chemkin_file(os.path.join(old_dir, 'chem_annotated.inp'), - os.path.join(old_dir, 'species_dictionary.txt'), - old_transport_path) - - new_species_list, new_reaction_list = load_chemkin_file(os.path.join(new_dir, 'chem_annotated.inp'), - os.path.join(new_dir, 'species_dictionary.txt'), - new_transport_path) + # define chemkin file path for old and new models + old_chemkin_path = os.path.join(old_dir, 'chem_annotated.inp') + new_chemkin_path = os.path.join(new_dir, 'chem_annotated.inp') + old_species_dict_path = os.path.join(old_dir, 'species_dictionary.txt') + new_species_dict_path = os.path.join(new_dir, 'species_dictionary.txt') + + surface = False + if os.path.exists(os.path.join(old_dir, 'chem_annotated-surface.inp')) and \ + os.path.exists(os.path.join(old_dir, 'chem_annotated-gas.inp')): + surface = True + old_chemkin_path = os.path.join(old_dir, 'chem_annotated-gas.inp') + new_chemkin_path = os.path.join(new_dir, 'chem_annotated-gas.inp') + old_surface_chemkin_path = os.path.join(old_dir, 'chem_annotated-surface.inp') + new_surface_chemkin_path = os.path.join(new_dir, 'chem_annotated-surface.inp') + ck2cti = True + + old_species_list, old_reaction_list = load_chemkin_file( + old_chemkin_path, + old_species_dict_path, + old_transport_path + ) + new_species_list, new_reaction_list = load_chemkin_file( + new_chemkin_path, + new_species_dict_path, + new_transport_path + ) + + old_surface_species_list = None + new_surface_species_list = None + new_surface_reaction_list = None + old_surface_reaction_list = None + if surface: + old_surface_species_list, old_surface_reaction_list = load_chemkin_file( + old_surface_chemkin_path, + old_species_dict_path, + old_transport_path + ) + new_surface_species_list, new_surface_reaction_list = load_chemkin_file( + new_surface_chemkin_path, + new_species_dict_path, + new_transport_path + ) self.old_sim = Cantera(species_list=old_species_list, - reaction_list=old_reaction_list, - output_directory=old_dir) + reaction_list=old_reaction_list, + output_directory=old_dir, + surface_species_list=old_surface_species_list, + surface_reaction_list=old_surface_reaction_list, + ) self.new_sim = Cantera(species_list=new_species_list, - reaction_list=new_reaction_list, - output_directory=new_dir) + reaction_list=new_reaction_list, + output_directory=new_dir, + surface_species_list=new_surface_species_list, + surface_reaction_list=new_surface_reaction_list, + ) # load each chemkin file into the cantera model if not ck2cti: self.old_sim.load_model() self.new_sim.load_model() else: - self.old_sim.load_chemkin_model(os.path.join(old_dir, 'chem_annotated.inp'), + self.old_sim.load_chemkin_model(old_chemkin_path, transport_file=old_transport_path, + surface_file=old_surface_chemkin_path, quiet=True) - self.new_sim.load_chemkin_model(os.path.join(new_dir, 'chem_annotated.inp'), + self.new_sim.load_chemkin_model(new_chemkin_path, transport_file=new_transport_path, + surface_file=new_surface_chemkin_path, quiet=True) def __str__(self): @@ -146,7 +188,7 @@ def __str__(self): """ return 'Observables Test Case: {0}'.format(self.title) - def generate_conditions(self, reactor_type_list, reaction_time_list, mol_frac_list, + def generate_conditions(self, reactor_type_list, reaction_time_list, mol_frac_list, surface_mol_frac_list=None, Tlist=None, Plist=None, Vlist=None): """ Creates a list of conditions from from the lists provided. @@ -161,6 +203,7 @@ def generate_conditions(self, reactor_type_list, reaction_time_list, mol_frac_li `reaction_time_list` A tuple object giving the ([list of reaction times], units) `mol_frac_list` A list of molfrac dictionaries with species object keys and mole fraction values + `surface_mol_frac_list` A list of molfrac dictionaries with species object keys for the surface To specify the system for an ideal gas, you must define 2 of the following 3 parameters: `T0List` A tuple giving the ([list of initial temperatures], units) 'P0List' A tuple giving the ([list of initial pressures], units) @@ -170,8 +213,12 @@ def generate_conditions(self, reactor_type_list, reaction_time_list, mol_frac_li """ # Store the conditions in the observables test case, for bookkeeping self.conditions = generate_cantera_conditions(reactor_type_list, reaction_time_list, mol_frac_list, + surface_mol_frac_list=surface_mol_frac_list, Tlist=Tlist, Plist=Plist, Vlist=Vlist) + if surface_mol_frac_list is None: + surface_mol_frac_list = [] # initialize here as list to avoid mutable default argument in function definition + # Map the mole fractions dictionaries to species objects from the old and new models old_mol_frac_list = [] new_mol_frac_list = [] @@ -192,12 +239,37 @@ def generate_conditions(self, reactor_type_list, reaction_time_list, mol_frac_li old_mol_frac_list.append(old_condition) new_mol_frac_list.append(new_condition) + if self.old_sim.surface: + old_surf_mol_frac_list = [] + new_surf_mol_frac_list = [] + for surf_mol_frac in surface_mol_frac_list: + old_surface_condition = {} + new_surface_condition = {} + old_surface_species_dict = get_rmg_species_from_user_species(list(surf_mol_frac.keys()), self.old_sim.surface_species_list) + new_surface_species_dict = get_rmg_species_from_user_species(list(surf_mol_frac.keys()), self.new_sim.surface_species_list) + for smiles, molfrac in surf_mol_frac.items(): + if old_surface_species_dict[smiles] is None: + raise Exception('SMILES {0} was not found in the old model!'.format(smiles)) + if new_surface_species_dict[smiles] is None: + raise Exception('SMILES {0} was not found in the new model!'.format(smiles)) + + old_surface_condition[old_surface_species_dict[smiles]] = molfrac + new_surface_condition[new_surface_species_dict[smiles]] = molfrac + old_surf_mol_frac_list.append(old_surface_condition) + new_surf_mol_frac_list.append(new_surface_condition) + else: + old_surf_mol_frac_list = [None] + new_surf_mol_frac_list = [None] + # Generate the conditions in each simulation self.old_sim.generate_conditions(reactor_type_list, reaction_time_list, old_mol_frac_list, + surface_mol_frac_list=old_surf_mol_frac_list, Tlist=Tlist, Plist=Plist, Vlist=Vlist) self.new_sim.generate_conditions(reactor_type_list, reaction_time_list, new_mol_frac_list, + surface_mol_frac_list=new_surf_mol_frac_list, Tlist=Tlist, Plist=Plist, Vlist=Vlist) + def compare(self, tol, plot=False): """ Compare the old and new model @@ -205,7 +277,7 @@ def compare(self, tol, plot=False): `plot`: if set to True, it will comparison plots of the two models comparing their species. Returns a list of variables failed in a list of tuples in the format: - + (CanteraCondition, variable label, variable_old, variable_new) """ @@ -220,10 +292,14 @@ def compare(self, tol, plot=False): print('{0} Comparison'.format(self)) # Check the species profile observables if 'species' in self.observables: - old_species_dict = get_rmg_species_from_user_species(self.observables['species'], self.old_sim.species_list) - new_species_dict = get_rmg_species_from_user_species(self.observables['species'], self.new_sim.species_list) - - # Check state variable observables + if self.old_sim.surface_species_list: + old_species_dict = get_rmg_species_from_user_species(self.observables['species'], self.old_sim.species_list + self.old_sim.surface_species_list) + new_species_dict = get_rmg_species_from_user_species(self.observables['species'], self.new_sim.species_list + self.new_sim.surface_species_list) + else: + old_species_dict = get_rmg_species_from_user_species(self.observables['species'], self.old_sim.species_list) + new_species_dict = get_rmg_species_from_user_species(self.observables['species'], self.new_sim.species_list) + + # Check state variable observables implemented_variables = ['temperature', 'pressure'] if 'variable' in self.observables: for item in self.observables['variable']: diff --git a/rmgpy/tools/regression.py b/rmgpy/tools/regression.py index a7c6e710ca..428918ae74 100644 --- a/rmgpy/tools/regression.py +++ b/rmgpy/tools/regression.py @@ -101,6 +101,15 @@ def read_input_file(path): new_imfs.append(new_dict) setups[3] = new_imfs + # convert keys of the initial surface mole fraction dictionaries (ismfs) from labels into Species objects. + if setups[5][0]: + ismfs = setups[5] + new_ismfs = [] + for ismf in ismfs: + new_dict = convert(ismf, initialSpecies) + new_ismfs.append(new_dict) + setups[5] = new_ismfs + return casetitle, observables, setups, tol @@ -115,14 +124,17 @@ def species(label, structure): return spc -def reactorSetups(reactorTypes, temperatures, pressures, initialMoleFractionsList, terminationTimes): +def reactorSetups(reactorTypes, temperatures, pressures, initialMoleFractionsList, terminationTimes, initialSurfaceMoleFractionsList=None): global setups + if initialSurfaceMoleFractionsList is None: + initialSurfaceMoleFractionsList = [None] + terminationTimes = Quantity(terminationTimes) temperatures = Quantity(temperatures) pressures = Quantity(pressures) - setups = [reactorTypes, temperatures, pressures, initialMoleFractionsList, terminationTimes] + setups = [reactorTypes, temperatures, pressures, initialMoleFractionsList, terminationTimes, initialSurfaceMoleFractionsList] def options(title='', tolerance=0.05): @@ -154,11 +166,12 @@ def run(benchmarkDir, testDir, title, observables, setups, tol): ck2cti=False, ) - reactor_types, temperatures, pressures, initial_mole_fractions_list, termination_times = setups + reactor_types, temperatures, pressures, initial_mole_fractions_list, termination_times, initial_surface_mole_fractions_list = setups case.generate_conditions( reactor_type_list=reactor_types, reaction_time_list=termination_times, mol_frac_list=initial_mole_fractions_list, + surface_mol_frac_list=initial_surface_mole_fractions_list, Tlist=temperatures, Plist=pressures ) From dad5e2790088ab36e1d4ac032b59f4b2835fe1af Mon Sep 17 00:00:00 2001 From: Sevy Harris Date: Mon, 6 Nov 2023 13:40:59 -0500 Subject: [PATCH 2/6] update CI workflow to include surface regression test --- .github/workflows/CI.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index ded6d65b51..75deef0616 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -167,7 +167,7 @@ jobs: id: regression-execution timeout-minutes: 60 run: | - for regr_test in aromatics liquid_oxidation nitrogen oxidation sulfur superminimal RMS_constantVIdealGasReactor_superminimal RMS_CSTR_liquid_oxidation RMS_liquidSurface_ch4o2cat fragment RMS_constantVIdealGasReactor_fragment; + for regr_test in aromatics liquid_oxidation nitrogen oxidation sulfur superminimal RMS_constantVIdealGasReactor_superminimal RMS_CSTR_liquid_oxidation RMS_liquidSurface_ch4o2cat fragment RMS_constantVIdealGasReactor_fragment minimal_surface; do if python-jl rmg.py test/regression/"$regr_test"/input.py; then echo "$regr_test" "Executed Successfully" @@ -242,7 +242,7 @@ jobs: run: | exec 2> >(tee -a regression.stderr >&2) 1> >(tee -a regression.stdout) mkdir -p "test/regression-diff" - for regr_test in aromatics liquid_oxidation nitrogen oxidation sulfur superminimal RMS_constantVIdealGasReactor_superminimal RMS_CSTR_liquid_oxidation fragment RMS_constantVIdealGasReactor_fragment; + for regr_test in aromatics liquid_oxidation nitrogen oxidation sulfur superminimal RMS_constantVIdealGasReactor_superminimal RMS_CSTR_liquid_oxidation fragment RMS_constantVIdealGasReactor_fragment minimal_surface; do echo "" echo "### Regression test $regr_test:" From e990de7386c250bce7e6c5a013d3885ab81068bf Mon Sep 17 00:00:00 2001 From: Sevy Harris Date: Wed, 2 Aug 2023 11:29:26 -0400 Subject: [PATCH 3/6] add input files for surface combustion test and example --- examples/rmg/minimal_surface/input.py | 146 ++++++++++++------ test/regression/minimal_surface/input.py | 123 +++++++++++++++ .../minimal_surface/regression_input.py | 67 ++++++++ 3 files changed, 291 insertions(+), 45 deletions(-) create mode 100644 test/regression/minimal_surface/input.py create mode 100644 test/regression/minimal_surface/regression_input.py diff --git a/examples/rmg/minimal_surface/input.py b/examples/rmg/minimal_surface/input.py index 5b45e1ee3d..7ef67e9502 100644 --- a/examples/rmg/minimal_surface/input.py +++ b/examples/rmg/minimal_surface/input.py @@ -1,67 +1,123 @@ # Data sources database( - thermoLibraries = ['primaryThermoLibrary'], - reactionLibraries = [], + thermoLibraries=['surfaceThermoPt111', 'primaryThermoLibrary', 'thermo_DFT_CCSDTF12_BAC','DFT_QCI_thermo'], # 'surfaceThermoPt' is the default. Thermo data is derived using bindingEnergies for other metals + reactionLibraries = [('Surface/CPOX_Pt/Deutschmann2006_adjusted', False)], # when Ni is used change the library to Surface/Deutschmann_Ni seedMechanisms = [], kineticsDepositories = ['training'], - kineticsFamilies = 'default', + kineticsFamilies = ['surface', 'default'], kineticsEstimator = 'rate rules', ) - -# List of species +catalystProperties( + bindingEnergies = { + 'H': (-2.75368,'eV/molecule'), + 'C': (-7.02516,'eV/molecule'), + 'N': (-4.63225,'eV/molecule'), + 'O': (-3.81153,'eV/molecule'), + }, + surfaceSiteDensity=(2.483e-9, 'mol/cm^2'), + coverageDependence=False +) species( - label='ethane', + label='CH4', reactive=True, - structure=SMILES("CC"), + structure=SMILES("C"), ) - -# Reaction systems -simpleReactor( - temperature=(1350,'K'), - pressure=(1.0,'bar'), - initialMoleFractions={ - "ethane": 1.0, - }, - terminationConversion={ - 'ethane': 0.9, - }, - terminationTime=(1e6,'s'), +species( + label='O2', + reactive=True, + structure=adjacencyList( + """ +1 O u1 p2 c0 {2,S} +2 O u1 p2 c0 {1,S} +"""), ) - -simpleReactor( - temperature=(1000,'K'), - pressure=(1.0,'bar'), - initialMoleFractions={ - "ethane": 1.0, +species( + label='N2', + reactive=False, + structure=SMILES("N#N"), +) +species( + label='X', + reactive=True, + structure=adjacencyList("1 X u0"), +) +# added for training +species( + label='CO2X', + reactive=True, + structure=adjacencyList( + """ + 1 O u0 p2 c0 {3,D} + 2 O u0 p2 c0 {3,D} + 3 C u0 p0 c0 {1,D} {2,D} + 4 X u0 p0 c0 + """), +) +species( + label='COX', + reactive=True, + structure=adjacencyList( + """ + 1 O u0 p2 c0 {2,D} + 2 C u0 p0 c0 {1,D} {3,D} + 3 X u0 p0 c0 {2,D} + """), +) +species( + label='OX', + reactive=True, + structure=adjacencyList( + """ + 1 O u0 p2 c0 {2,D} + 2 X u0 p0 c0 {1,D} + """), +) +# If you would like to forbid the bidentate form of absorbed CO2 from your model, +# use the following `CO2_bidentate` forbidden structure +# forbidden( +# label='CO2_bidentate', +# structure=SMILES("O=C(*)O*"), +# ) +#---------- +# Reaction systems +surfaceReactor( + temperature=(1000, 'K'), + initialPressure=(1.0, 'bar'), + initialGasMoleFractions={ + "CH4": 0.15, + "O2": 0.15, + "N2": 0.7, }, - terminationConversion={ - 'ethane': 0.9, + initialSurfaceCoverages={ + "X": 1.0, }, - terminationTime=(1e6,'s'), + surfaceVolumeRatio=(1.e5, 'm^-1'), + terminationConversion = { "CH4": 0.95,}, + terminationTime=(0.1, 's'), + terminationRateRatio=0.01, ) - - - simulator( - atol=1e-16, - rtol=1e-8, + atol=1e-18, + rtol=1e-12, ) - model( toleranceKeepInEdge=0.0, - toleranceMoveToCore=100.0, - toleranceInterruptSimulation=100.0, + toleranceMoveToCore=1e-1, + toleranceInterruptSimulation=0.1, maximumEdgeSpecies=100000, - toleranceMoveEdgeReactionToSurface=1.0, - toleranceMoveSurfaceReactionToCore=2.0, - toleranceMoveSurfaceSpeciesToCore=.001, - dynamicsTimeScale=(1.0e-10,'sec'), ) - options( units='si', - generateOutputHTML=True, - generatePlots=False, - saveEdgeSpecies=True, - saveSimulationProfiles=True, + generateOutputHTML=False, + generatePlots=False, # Enable to make plots of core and edge size etc. But takes a lot of the total runtime! + saveEdgeSpecies=False, + saveSimulationProfiles=False, +) +generatedSpeciesConstraints( + allowed=['input species','reaction libraries'], + maximumCarbonAtoms=2, + maximumOxygenAtoms=2, + maximumSurfaceSites=2, ) + + diff --git a/test/regression/minimal_surface/input.py b/test/regression/minimal_surface/input.py new file mode 100644 index 0000000000..a2d8c101de --- /dev/null +++ b/test/regression/minimal_surface/input.py @@ -0,0 +1,123 @@ +# Data sources +database( + thermoLibraries=['surfaceThermoPt111', 'primaryThermoLibrary', 'thermo_DFT_CCSDTF12_BAC','DFT_QCI_thermo'], # 'surfaceThermoPt' is the default. Thermo data is derived using bindingEnergies for other metals + reactionLibraries = [('Surface/CPOX_Pt/Deutschmann2006_adjusted', False)], # when Ni is used change the library to Surface/Deutschmann_Ni + seedMechanisms = [], + kineticsDepositories = ['training'], + kineticsFamilies = ['surface', 'default'], + kineticsEstimator = 'rate rules', +) +catalystProperties( + bindingEnergies = { + 'H': (-2.75368,'eV/molecule'), + 'C': (-7.02516,'eV/molecule'), + 'N': (-4.63225,'eV/molecule'), + 'O': (-3.81153,'eV/molecule'), + }, + surfaceSiteDensity=(2.483e-9, 'mol/cm^2'), + coverageDependence=False +) +species( + label='CH4', + reactive=True, + structure=SMILES("C"), +) +species( + label='O2', + reactive=True, + structure=adjacencyList( + """ +1 O u1 p2 c0 {2,S} +2 O u1 p2 c0 {1,S} +"""), +) +species( + label='N2', + reactive=False, + structure=SMILES("N#N"), +) +species( + label='X', + reactive=True, + structure=adjacencyList("1 X u0"), +) +# added for training +species( + label='CO2X', + reactive=True, + structure=adjacencyList( + """ + 1 O u0 p2 c0 {3,D} + 2 O u0 p2 c0 {3,D} + 3 C u0 p0 c0 {1,D} {2,D} + 4 X u0 p0 c0 + """), +) +species( + label='COX', + reactive=True, + structure=adjacencyList( + """ + 1 O u0 p2 c0 {2,D} + 2 C u0 p0 c0 {1,D} {3,D} + 3 X u0 p0 c0 {2,D} + """), +) +species( + label='OX', + reactive=True, + structure=adjacencyList( + """ + 1 O u0 p2 c0 {2,D} + 2 X u0 p0 c0 {1,D} + """), +) +# If you would like to forbid the bidentate form of absorbed CO2 from your model, +# use the following `CO2_bidentate` forbidden structure +# forbidden( +# label='CO2_bidentate', +# structure=SMILES("O=C(*)O*"), +# ) +#---------- +# Reaction systems +surfaceReactor( + temperature=(1000, 'K'), + initialPressure=(1.0, 'bar'), + initialGasMoleFractions={ + "CH4": 0.15, + "O2": 0.15, + "N2": 0.7, + }, + initialSurfaceCoverages={ + "X": 1.0, + }, + surfaceVolumeRatio=(1.e5, 'm^-1'), + terminationConversion={"CH4": 0.95,}, + terminationTime=(0.1, 's'), + terminationRateRatio=0.01, +) +simulator( + atol=1e-18, + rtol=1e-12, +) +model( + toleranceKeepInEdge=0.0, + toleranceMoveToCore=1e-1, + toleranceInterruptSimulation=0.1, + maximumEdgeSpecies=100000, +) +options( + units='si', + generateOutputHTML=False, + generatePlots=False, # Enable to make plots of core and edge size etc. But takes a lot of the total runtime! + saveEdgeSpecies=False, + saveSimulationProfiles=False, +) +generatedSpeciesConstraints( + allowed=['input species','reaction libraries'], + maximumCarbonAtoms=2, + maximumOxygenAtoms=2, + maximumSurfaceSites=2, +) + + diff --git a/test/regression/minimal_surface/regression_input.py b/test/regression/minimal_surface/regression_input.py new file mode 100644 index 0000000000..785705e38c --- /dev/null +++ b/test/regression/minimal_surface/regression_input.py @@ -0,0 +1,67 @@ +options( + title = 'minimal_surface', + tolerance = 0.5, +) + +# observables +observable( + label = 'CH4', + structure=SMILES('C'), +) + +observable( + label = 'O2', + structure=SMILES('[O][O]'), +) + +observable( + label = 'X', + structure=adjacencyList( + """ +1 X u0 p0 c0 +"""), +) + + +# List of species +species( + label='CH4', + structure=SMILES("[CH4]"), +) +species( + label='O2', + structure=adjacencyList( + """ +1 O u1 p2 c0 {2,S} +2 O u1 p2 c0 {1,S} +"""), +) +species( + label='N2', + structure=SMILES("N#N"), +) +species( + label='X', + structure=adjacencyList( + """ +1 X u0 p0 c0 +"""), +) + + +# reactor setups +reactorSetups( + reactorTypes=['IdealGasReactor'], + terminationTimes=([1.0], 's'), + initialMoleFractionsList=[{ + 'CH4': 0.15, + 'O2': 0.15, + 'N2': 0.7, + }], + initialSurfaceMoleFractionsList=[{ + 'X': 1.0, + }], + + temperatures=([1000.0], 'K'), + pressures=([1.0], 'bar'), +) From 6a149da762e166cd105cc0f2484411585a957981 Mon Sep 17 00:00:00 2001 From: Sevy Harris Date: Fri, 3 May 2024 14:36:57 -0400 Subject: [PATCH 4/6] add workaround to check for surface files if calling function assumed gas-phase --- rmgpy/tools/diffmodels.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/rmgpy/tools/diffmodels.py b/rmgpy/tools/diffmodels.py index 51c927493e..e93609a77e 100644 --- a/rmgpy/tools/diffmodels.py +++ b/rmgpy/tools/diffmodels.py @@ -344,6 +344,23 @@ def execute(chemkin1, species_dict1, thermo1, chemkin2, species_dict2, thermo2, model1 = ReactionModel() model2 = ReactionModel() + # check if this is supposed to be comparing surface mechanisms + chemkin_gas1 = chemkin1.replace('.inp', '-gas.inp') + chemkin_surface1 = chemkin1.replace('.inp', '-surface.inp') + if not os.path.exists(chemkin1) and \ + os.path.exists(chemkin_surface1) and \ + os.path.exists(chemkin_gas1): + chemkin1 = chemkin_gas1 + kwargs['surface_path1'] = chemkin_surface1 + + chemkin_gas2 = chemkin2.replace('.inp', '-gas.inp') + chemkin_surface2 = chemkin2.replace('.inp', '-surface.inp') + if not os.path.exists(chemkin2) and \ + os.path.exists(chemkin_surface2) and \ + os.path.exists(chemkin_gas2): + chemkin2 = chemkin_gas2 + kwargs['surface_path2'] = chemkin_surface2 + try: surface_path1 = kwargs['surface_path1'] surface_path2 = kwargs['surface_path2'] From afd00ae0281b5b120caa8f11f4d29bdc03111492 Mon Sep 17 00:00:00 2001 From: Sevy Harris Date: Sun, 5 May 2024 16:15:40 -0400 Subject: [PATCH 5/6] turn edge species on for regression comparison --- test/regression/minimal_surface/input.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/regression/minimal_surface/input.py b/test/regression/minimal_surface/input.py index a2d8c101de..ddfea2dd87 100644 --- a/test/regression/minimal_surface/input.py +++ b/test/regression/minimal_surface/input.py @@ -110,7 +110,7 @@ units='si', generateOutputHTML=False, generatePlots=False, # Enable to make plots of core and edge size etc. But takes a lot of the total runtime! - saveEdgeSpecies=False, + saveEdgeSpecies=True, saveSimulationProfiles=False, ) generatedSpeciesConstraints( From 4cd90fb14961cc2444dda713dddba1dde9ef8c7a Mon Sep 17 00:00:00 2001 From: Sevy Harris Date: Mon, 6 May 2024 09:05:38 -0400 Subject: [PATCH 6/6] add a check for all zeros in thermo comparison --- rmgpy/thermo/model.pyx | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/rmgpy/thermo/model.pyx b/rmgpy/thermo/model.pyx index 7a4f0a663c..502c863925 100644 --- a/rmgpy/thermo/model.pyx +++ b/rmgpy/thermo/model.pyx @@ -163,13 +163,18 @@ cdef class HeatCapacityModel(RMGObject): Tdata = [300,400,500,600,800,1000,1500,2000] for T in Tdata: - if not (0.8 < self.get_heat_capacity(T) / other.get_heat_capacity(T) < 1.25): + # Do exact comparison in addition to relative in case both are zero (surface site) + if self.get_heat_capacity(T) != other.get_heat_capacity(T) and \ + not (0.8 < self.get_heat_capacity(T) / other.get_heat_capacity(T) < 1.25): return False - elif not (0.8 < self.get_enthalpy(T) / other.get_enthalpy(T) < 1.25): + elif self.get_enthalpy(T) != other.get_enthalpy(T) and \ + not (0.8 < self.get_enthalpy(T) / other.get_enthalpy(T) < 1.25): return False - elif not (0.8 < self.get_entropy(T) / other.get_entropy(T) < 1.25): + elif self.get_entropy(T) != other.get_entropy(T) and \ + not (0.8 < self.get_entropy(T) / other.get_entropy(T) < 1.25): return False - elif not (0.8 < self.get_free_energy(T) / other.get_free_energy(T) < 1.25): + elif self.get_free_energy(T) != other.get_free_energy(T) and \ + not (0.8 < self.get_free_energy(T) / other.get_free_energy(T) < 1.25): return False return True @@ -185,13 +190,18 @@ cdef class HeatCapacityModel(RMGObject): Tdata = [300,400,500,600,800,1000,1500,2000] for T in Tdata: - if not (0.95 < self.get_heat_capacity(T) / other.get_heat_capacity(T) < 1.05): + # Do exact comparison in addition to relative in case both are zero (surface site) + if self.get_heat_capacity(T) != other.get_heat_capacity(T) and \ + not (0.95 < self.get_heat_capacity(T) / other.get_heat_capacity(T) < 1.05): return False - elif not (0.95 < self.get_enthalpy(T) / other.get_enthalpy(T) < 1.05): + elif self.get_enthalpy(T) != other.get_enthalpy(T) and \ + not (0.95 < self.get_enthalpy(T) / other.get_enthalpy(T) < 1.05): return False - elif not (0.95 < self.get_entropy(T) / other.get_entropy(T) < 1.05): + elif self.get_entropy(T) != other.get_entropy(T) and \ + not (0.95 < self.get_entropy(T) / other.get_entropy(T) < 1.05): return False - elif not (0.95 < self.get_free_energy(T) / other.get_free_energy(T) < 1.05): + elif self.get_free_energy(T) != other.get_free_energy(T) and \ + not (0.95 < self.get_free_energy(T) / other.get_free_energy(T) < 1.05): return False return True