From 67f934da2aebf8fa07d1d1a5049270fdff9e1316 Mon Sep 17 00:00:00 2001 From: Sevy Harris Date: Thu, 19 Oct 2023 14:44:59 -0400 Subject: [PATCH 01/22] 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 02/22] 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 03/22] 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 04/22] 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 05/22] 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 06/22] 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 From 55bfc6558638d1dbc99ca69b7ab5f5086e459aea Mon Sep 17 00:00:00 2001 From: Richard West Date: Wed, 15 May 2024 14:39:57 -0400 Subject: [PATCH 07/22] Better debugging messages for database errors. When there are issues in a database file, such as nan values here: uncertainty=RateUncertainty(mu=nan, var=nan, which may be due to a bug in the tree fitting, it was very hard to find the cause of the error because the stack trace wouldn't tell you where the error was, and the file handle had been read by the time you caught the exception. With this change, we read the file to a variable, then try executing it, and catch and report any exceptions properly. --- rmgpy/data/base.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/rmgpy/data/base.py b/rmgpy/data/base.py index 872a1d18d6..54c0e6c070 100644 --- a/rmgpy/data/base.py +++ b/rmgpy/data/base.py @@ -236,13 +236,17 @@ def load(self, path, local_context=None, global_context=None): local_context[key] = value # Process the file - f = open(path, 'r') + with open(path, 'r') as f: + content = f.read() try: - exec(f.read(), global_context, local_context) - except Exception: - logging.error('Error while reading database {0!r}.'.format(path)) + exec(content, global_context, local_context) + except Exception as e: + logging.exception(f'Error while reading database file {path}.') + line_number = e.__traceback__.tb_next.tb_lineno + logging.error(f'Error occurred at or near line {line_number} of {path}.') + lines = content.splitlines() + logging.error(f'Line: {lines[line_number - 1]}') raise - f.close() # Extract the database metadata self.name = local_context['name'] From 192daeaeb35dfb155a30969f4da22899f4c4bc80 Mon Sep 17 00:00:00 2001 From: Jackson Burns Date: Fri, 7 Jun 2024 09:11:28 -0400 Subject: [PATCH 08/22] github actions CI maintainence - increment miniconda setup to v3 from v2 - set macos build runner to use intel macs rather than m1, since we don't have binaries for m1 - switch download artifact step to use the official github one, which now allows grabbing artifacts from other runs Resolves https://github.com/ReactionMechanismGenerator/RMG-Py/issues/2611 See https://github.com/ReactionMechanismGenerator/RMG-Py/pull/2667 for more history asdf --- .github/workflows/CI.yml | 33 +++++++++++++++++---------------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 75deef0616..71cdc6f287 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -51,7 +51,7 @@ env: jobs: build-osx: - runs-on: macos-latest + runs-on: macos-13 # skip scheduled runs from forks if: ${{ !( github.repository != 'ReactionMechanismGenerator/RMG-Py' && github.event_name == 'schedule' ) }} defaults: @@ -63,7 +63,7 @@ jobs: # configures the mamba environment manager and builds the environment - name: Setup Mambaforge Python 3.7 - uses: conda-incubator/setup-miniconda@v2 + uses: conda-incubator/setup-miniconda@v3 with: environment-file: environment.yml miniforge-variant: Mambaforge @@ -113,7 +113,7 @@ jobs: # configures the mamba environment manager and builds the environment - name: Setup Mambaforge Python 3.7 - uses: conda-incubator/setup-miniconda@v2 + uses: conda-incubator/setup-miniconda@v3 with: environment-file: environment.yml miniforge-variant: Mambaforge @@ -185,7 +185,7 @@ jobs: # Upload Regression Results as Failed if above step failed - name: Upload Failed Results if: ${{ failure() && steps.regression-execution.conclusion == 'failure' }} - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@v4 with: name: failed_regression_results path: | @@ -195,7 +195,7 @@ jobs: - name: Upload Results as Reference # upload the results for scheduled CI (on main) and pushes to main if: ${{ env.REFERENCE_JOB == 'true' }} - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@v4 with: name: stable_regression_results path: | @@ -204,7 +204,7 @@ jobs: # Upload Regression Results as Dynamic if Push to non-main Branch - name: Upload Results as Dynamic if: ${{ env.REFERENCE_JOB == 'false' }} - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@v4 with: name: dynamic_regression_results path: | @@ -215,23 +215,24 @@ jobs: run: mkdir stable_regression_results # Retrieve Stable Results for reference - # Will need to use this -> https://github.com/dawidd6/action-download-artifact + - name : Find ID of Reference Results + env: + GH_TOKEN: ${{ secrets.GITHUB_TOKEN }} + run: | + run_id=$(gh run list -R ReactionMechanismGenerator/RMG-Py --workflow="Continuous Integration" --branch main --limit 1 --json databaseId --jq '.[0].databaseId') + echo "CI_RUN_ID=$run_id" >> $GITHUB_ENV + - name: Retrieve Stable Regression Results if: ${{ env.REFERENCE_JOB == 'false' }} - uses: dsnopek/action-download-artifact@91dda23aa09c68860977dd0ed11d93c0ed3795e7 # see https://github.com/ReactionMechanismGenerator/RMG-Py/pull/2459#issuecomment-1582850815 + uses: actions/download-artifact@v4 with: # this will search for the last successful execution of CI on main and download # the stable regression results - workflow: CI.yml - workflow_conclusion: success - repo: ReactionMechanismGenerator/RMG-Py - branch: main + run-id: ${{ env.CI_RUN_ID }} + repository: ReactionMechanismGenerator/RMG-Py + github-token: ${{ secrets.GH_PAT }} name: stable_regression_results path: stable_regression_results - search_artifacts: true # retrieves the last run result, either scheduled daily or on push to main - ensure_latest: true # ensures that the latest run is retrieved - # should result in a set of folders inside stable_regression_results - # each of which has the stable result for that example/test # Regression Testing - Actual Comparisons - name: Regression Tests - Compare to Baseline From 6cba2b74c1018eccaa87cbb5bdff4e6c0aa7aaeb Mon Sep 17 00:00:00 2001 From: donerancl Date: Wed, 31 Jan 2024 17:05:12 -0500 Subject: [PATCH 09/22] squash --- arkane/ess/gaussian.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/arkane/ess/gaussian.py b/arkane/ess/gaussian.py index 2e84d6f231..f00d9dd9e0 100644 --- a/arkane/ess/gaussian.py +++ b/arkane/ess/gaussian.py @@ -35,7 +35,7 @@ import logging import math import os.path - +import re import numpy as np import rmgpy.constants as constants @@ -309,9 +309,8 @@ def load_energy(self, zpe_scale_factor=1.): with open(self.path, 'r') as f: line = f.readline() while line != '': - if 'SCF Done:' in line: - e_elect = float(line.split()[4]) * constants.E_h * constants.Na + e_elect = float(re.findall(r"SCF Done: E\(.+\) \=\s+[^\s]+",line)[0].split()[-1])* constants.E_h * constants.Na elect_energy_source = 'SCF' elif ' E2(' in line and ' E(' in line: e_elect = float(line.split()[-1].replace('D', 'E')) * constants.E_h * constants.Na @@ -351,7 +350,7 @@ def load_energy(self, zpe_scale_factor=1.): # G4MP2 calculation without opt and freq calculation # Keyword in Gaussian G4MP2(SP), No zero-point or thermal energies are included. e_elect = float(line.split()[2]) * constants.E_h * constants.Na - + # Read the ZPE from the "E(ZPE)=" line, as this is the scaled version. # Gaussian defines the following as # E (0 K) = Elec + E(ZPE), @@ -376,6 +375,12 @@ def load_energy(self, zpe_scale_factor=1.): elect_energy_source = 'HF' except ValueError: pass + elif 'Energy=' in line: + # for xtb + e_elect = float(line.split()[1]) * constants.E_h * constants.Na + elif 'Energy=' in line: + # for xtb + e_elect = float(line.split()[1]) * constants.E_h * constants.Na # Read the next line in the file line = f.readline() From 34e3d5d19c3febd76095e103e56708ce6421f8f8 Mon Sep 17 00:00:00 2001 From: Anna Doner <47831564+donerancl@users.noreply.github.com> Date: Mon, 5 Feb 2024 15:07:52 -0500 Subject: [PATCH 10/22] remove duplicate lines raising error I think I added these in accidentally, I got the same error and omitting them fixed it --- arkane/ess/gaussian.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/arkane/ess/gaussian.py b/arkane/ess/gaussian.py index f00d9dd9e0..da1b57ef2b 100644 --- a/arkane/ess/gaussian.py +++ b/arkane/ess/gaussian.py @@ -375,12 +375,6 @@ def load_energy(self, zpe_scale_factor=1.): elect_energy_source = 'HF' except ValueError: pass - elif 'Energy=' in line: - # for xtb - e_elect = float(line.split()[1]) * constants.E_h * constants.Na - elif 'Energy=' in line: - # for xtb - e_elect = float(line.split()[1]) * constants.E_h * constants.Na # Read the next line in the file line = f.readline() From a597078c8e11225b9b314d630d398ce59bbada4e Mon Sep 17 00:00:00 2001 From: Jackson Burns <33505528+JacksonBurns@users.noreply.github.com> Date: Fri, 7 Jun 2024 14:17:38 -0400 Subject: [PATCH 11/22] Apply suggestions from code review en-pretty-fy whitespace Co-authored-by: Hao-Wei Pang <45482070+hwpang@users.noreply.github.com> --- arkane/ess/gaussian.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/arkane/ess/gaussian.py b/arkane/ess/gaussian.py index da1b57ef2b..156823b987 100644 --- a/arkane/ess/gaussian.py +++ b/arkane/ess/gaussian.py @@ -310,7 +310,7 @@ def load_energy(self, zpe_scale_factor=1.): line = f.readline() while line != '': if 'SCF Done:' in line: - e_elect = float(re.findall(r"SCF Done: E\(.+\) \=\s+[^\s]+",line)[0].split()[-1])* constants.E_h * constants.Na + e_elect = float(re.findall(r"SCF Done: E\(.+\) \=\s+[^\s]+", line)[0].split()[-1]) * constants.E_h * constants.Na elect_energy_source = 'SCF' elif ' E2(' in line and ' E(' in line: e_elect = float(line.split()[-1].replace('D', 'E')) * constants.E_h * constants.Na @@ -350,7 +350,6 @@ def load_energy(self, zpe_scale_factor=1.): # G4MP2 calculation without opt and freq calculation # Keyword in Gaussian G4MP2(SP), No zero-point or thermal energies are included. e_elect = float(line.split()[2]) * constants.E_h * constants.Na - # Read the ZPE from the "E(ZPE)=" line, as this is the scaled version. # Gaussian defines the following as # E (0 K) = Elec + E(ZPE), From 25e04b5611459d04ee4aaf746fa827b15c24eda0 Mon Sep 17 00:00:00 2001 From: Jackson Burns Date: Sat, 8 Jun 2024 22:10:37 -0400 Subject: [PATCH 12/22] pass the token through with the correct variable name that's what i get for copy-pasting from the README --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 71cdc6f287..b9b8ed39ed 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -230,7 +230,7 @@ jobs: # the stable regression results run-id: ${{ env.CI_RUN_ID }} repository: ReactionMechanismGenerator/RMG-Py - github-token: ${{ secrets.GH_PAT }} + github-token: ${{ secrets.GITHUB_TOKEN }} name: stable_regression_results path: stable_regression_results From 4889ca09f6886ca3b4db73f508ffcc24df29c601 Mon Sep 17 00:00:00 2001 From: Jackson Burns <33505528+JacksonBurns@users.noreply.github.com> Date: Sun, 9 Jun 2024 06:48:01 -0400 Subject: [PATCH 13/22] temporarily remove `minimal_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 b9b8ed39ed..7b5a39c162 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 minimal_surface; + 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; do if python-jl rmg.py test/regression/"$regr_test"/input.py; then echo "$regr_test" "Executed Successfully" @@ -243,7 +243,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 minimal_surface; + for regr_test in aromatics liquid_oxidation nitrogen oxidation sulfur superminimal RMS_constantVIdealGasReactor_superminimal RMS_CSTR_liquid_oxidation fragment RMS_constantVIdealGasReactor_fragment; do echo "" echo "### Regression test $regr_test:" From b7761fd8b6fd746d0642d83287fe38c51ba17b9c Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Thu, 16 May 2024 21:17:48 +0300 Subject: [PATCH 14/22] Send P in SI units to Chebyshev.fit_to_data() when reversing a rate The Chebyshev.fit_to_data() method treats the Pmin and Pmax args as SI (Pa units), see here for example: `self.Pmin = (Pmin * 1e-5, "bar")`. Make sure we transmit the arguments in SI units and not in the arbitrary units they were defined in. --- rmgpy/reaction.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index 7484f70d43..f10e1dcefa 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -1030,15 +1030,15 @@ def generate_reverse_rate_coefficient(self, network_kinetics=False, Tmin=None, T return self.reverse_arrhenius_rate(kf, kunits) elif isinstance(kf, Chebyshev): - Tlist = 1.0 / np.linspace(1.0 / kf.Tmax.value, 1.0 / kf.Tmin.value, 50) - Plist = np.linspace(kf.Pmin.value, kf.Pmax.value, 20) + Tlist = 1.0 / np.linspace(1.0 / kf.Tmax.value_si, 1.0 / kf.Tmin.value_si, 50) + Plist = np.linspace(kf.Pmin.value_si, kf.Pmax.value_si, 20) K = np.zeros((len(Tlist), len(Plist)), float) for Tindex, T in enumerate(Tlist): for Pindex, P in enumerate(Plist): K[Tindex, Pindex] = kf.get_rate_coefficient(T, P) / self.get_equilibrium_constant(T) kr = Chebyshev() - kr.fit_to_data(Tlist, Plist, K, kunits, kf.degreeT, kf.degreeP, kf.Tmin.value, kf.Tmax.value, kf.Pmin.value, - kf.Pmax.value) + kr.fit_to_data(Tlist, Plist, K, kunits, kf.degreeT, kf.degreeP, kf.Tmin.value, kf.Tmax.value, + kf.Pmin.value_si, kf.Pmax.value_si) return kr elif isinstance(kf, PDepArrhenius): From f51edfafd12b97e2c126199b22e9663c135d0401 Mon Sep 17 00:00:00 2001 From: Jackson Burns Date: Fri, 14 Jun 2024 08:17:58 -0400 Subject: [PATCH 15/22] fix expected critical pressure see: https://github.com/ReactionMechanismGenerator/RMG-Py/pull/2682#issuecomment-2167058250 --- test/rmgpy/data/transportTest.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/rmgpy/data/transportTest.py b/test/rmgpy/data/transportTest.py index 9e3b9f71bf..0b97451142 100644 --- a/test/rmgpy/data/transportTest.py +++ b/test/rmgpy/data/transportTest.py @@ -153,7 +153,7 @@ def test_joback(self): "CC(=O)C", Length(5.36421, "angstroms"), Energy(3.20446, "kJ/mol"), - "Epsilon & sigma estimated with Tc=500.53 K, Pc=47.11 bar (from Joback method)", + "Epsilon & sigma estimated with Tc=500.53 K, Pc=48.02 bar (from Joback method)", ], [ "cyclopenta-1,2-diene", From f08312dbfb034ab9b9927e9614ac80fefdbaf30c Mon Sep 17 00:00:00 2001 From: Jackson Burns <33505528+JacksonBurns@users.noreply.github.com> Date: Fri, 14 Jun 2024 09:36:01 -0400 Subject: [PATCH 16/22] update expected sigma value --- test/rmgpy/data/transportTest.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/rmgpy/data/transportTest.py b/test/rmgpy/data/transportTest.py index 0b97451142..45d7855596 100644 --- a/test/rmgpy/data/transportTest.py +++ b/test/rmgpy/data/transportTest.py @@ -151,7 +151,7 @@ def test_joback(self): [ "acetone", "CC(=O)C", - Length(5.36421, "angstroms"), + Length(5.32979, "angstroms"), Energy(3.20446, "kJ/mol"), "Epsilon & sigma estimated with Tc=500.53 K, Pc=48.02 bar (from Joback method)", ], From e661630ec0a7b2a1cf0b5049612796b06f320476 Mon Sep 17 00:00:00 2001 From: Jackson Burns <33505528+JacksonBurns@users.noreply.github.com> Date: Sat, 15 Jun 2024 12:22:45 -0400 Subject: [PATCH 17/22] switch multiple `rmg` packages to their official versions --- environment.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/environment.yml b/environment.yml index 64df68a1c2..48a4da34a3 100644 --- a/environment.yml +++ b/environment.yml @@ -81,13 +81,13 @@ dependencies: - conda-forge::gprof2dot - conda-forge::numdifftools - conda-forge::quantities + - conda-forge::muq + - conda-forge::lpsolve55 + - conda-forge::ringdecomposerlib-python # packages we maintain - - rmg::lpsolve55 - - rmg::muq2 - rmg::pydas >=1.0.3 - rmg::pydqed >=1.0.3 - - rmg::pyrdl - rmg::pyrms - rmg::symmetry From c3c584d8b5c5baadc7e35ddbec542cc27b84c914 Mon Sep 17 00:00:00 2001 From: Jackson Burns Date: Sat, 15 Jun 2024 12:30:24 -0400 Subject: [PATCH 18/22] force make to use gcc and g++ compilers changed in Makefile to avoid clang issue on mac os, as clang is now installed with muq and may conflict during installation --- Makefile | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Makefile b/Makefile index c573fc3b47..a2f6b3fed2 100644 --- a/Makefile +++ b/Makefile @@ -4,6 +4,9 @@ # ################################################################################ +CC=gcc +CXX=g++ + .PHONY : all minimal main solver check pycheck arkane clean install decython documentation test q2dtor all: pycheck main solver check From efb8e69c1220ec41281120a33a8ee938ef848c72 Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Sat, 18 May 2024 10:15:28 +0300 Subject: [PATCH 19/22] Added external_library_labels attribute to KineticsDatabase --- rmgpy/data/kinetics/database.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/rmgpy/data/kinetics/database.py b/rmgpy/data/kinetics/database.py index 41af9a6b55..c3bf6b0804 100644 --- a/rmgpy/data/kinetics/database.py +++ b/rmgpy/data/kinetics/database.py @@ -62,6 +62,7 @@ def __init__(self): self.recommended_families = {} self.families = {} self.libraries = {} + self.external_library_labels = {} self.library_order = [] # a list of tuples in the format ('library_label', LibraryType), # where LibraryType is set to either 'Reaction Library' or 'Seed'. self.local_context = { @@ -227,7 +228,7 @@ def load_libraries(self, path, libraries=None): The `path` points to the folder of kinetics libraries in the database, and the libraries should be in files like :file:`/.py`. """ - + self.external_library_labels = dict() if libraries is not None: for library_name in libraries: library_file = os.path.join(path, library_name, 'reactions.py') @@ -238,6 +239,7 @@ def load_libraries(self, path, libraries=None): library = KineticsLibrary(label=short_library_name) library.load(library_file, self.local_context, self.global_context) self.libraries[library.label] = library + self.external_library_labels[library_name] = library.label elif os.path.exists(library_file): logging.info(f'Loading kinetics library {library_name} from {library_file}...') library = KineticsLibrary(label=library_name) From 97f19d0de08dff2a012879b57f082eb930cab5d6 Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Thu, 13 Jun 2024 08:10:08 +0300 Subject: [PATCH 20/22] Strip os.path.sep when extracting the name of an external library --- rmgpy/data/kinetics/database.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/data/kinetics/database.py b/rmgpy/data/kinetics/database.py index c3bf6b0804..91b13e3972 100644 --- a/rmgpy/data/kinetics/database.py +++ b/rmgpy/data/kinetics/database.py @@ -234,7 +234,7 @@ def load_libraries(self, path, libraries=None): library_file = os.path.join(path, library_name, 'reactions.py') if os.path.exists(library_name): library_file = os.path.join(library_name, 'reactions.py') - short_library_name = os.path.split(library_name)[-1] + short_library_name = os.path.basename(library_name.rstrip(os.path.sep)) logging.info(f'Loading kinetics library {short_library_name} from {library_name}...') library = KineticsLibrary(label=short_library_name) library.load(library_file, self.local_context, self.global_context) From 2b153a7fdc566c612975828bd5b32d310cc480da Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Thu, 13 Jun 2024 07:23:11 +0300 Subject: [PATCH 21/22] Replace an external reaction library path with the library label when loading a library, since kinetic libraries are stored in database.kinetics.libraries by their labels --- rmgpy/rmg/model.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index 1d8c581cc0..1ab45dc67c 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -1718,7 +1718,12 @@ def add_reaction_library_to_edge(self, reaction_library): num_old_edge_reactions = len(self.edge.reactions) logging.info("Adding reaction library {0} to model edge...".format(reaction_library)) - reaction_library = database.kinetics.libraries[reaction_library] + if reaction_library in database.kinetics.libraries: + reaction_library = database.kinetics.libraries[reaction_library] + elif reaction_library in database.kinetics.external_library_labels: + reaction_library = database.kinetics.libraries[database.kinetics.external_library_labels[reaction_library]] + else: + raise ValueError(f'Library {reaction_library} not found.') rxns = reaction_library.get_library_reactions() for rxn in rxns: From 6d9a56e636648b93e330fd0afaa294f23c148103 Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Sat, 25 May 2024 13:52:30 +0300 Subject: [PATCH 22/22] Docs: Documented the external library feature --- documentation/source/users/rmg/input.rst | 21 +++++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/documentation/source/users/rmg/input.rst b/documentation/source/users/rmg/input.rst index 099a824c11..38dead4517 100644 --- a/documentation/source/users/rmg/input.rst +++ b/documentation/source/users/rmg/input.rst @@ -49,7 +49,7 @@ by Benson's method. For example, if you wish to use the GRI-Mech 3.0 mechanism [GRIMech3.0]_ as a ThermoLibrary in your model, the syntax will be:: - thermoLibraries = ['primaryThermoLibrary','GRI-Mech3.0'] + thermoLibraries = ['primaryThermoLibrary', 'GRI-Mech3.0'] .. [GRIMech3.0] Gregory P. Smith, David M. Golden, Michael Frenklach, Nigel W. Moriarty, Boris Eiteneer, Mikhail Goldenberg, C. Thomas Bowman, Ronald K. Hanson, Soonho Song, William C. Gardiner, Jr., Vitali V. Lissianski, and Zhiwei Qin http://combustion.berkeley.edu/gri-mech/ @@ -78,7 +78,7 @@ In the following example, the user has created a reaction library with a few additional reactions specific to n-butane, and these reactions are to be used in addition to the Glarborg C3 library:: - reactionLibraries = [('Glarborg/C3',False)], + reactionLibraries = [('Glarborg/C3', False)], The keyword False/True permits user to append all unused reactions (= kept in the edge) from this library to the chemkin file. True means those reactions will be appended. Using just the string inputs would lead to @@ -104,6 +104,23 @@ given in each mechanism, the different mechanisms can have different units. Library. +.. _externallib: + +External Libraries +------------------ +Users may direct RMG to use thermo and/or kinetic libraries which are not included in the RMG database, +e.g., a library a user created that was intentionally saved to a path different than the conventional +RMG-database location. In such cases, the user can specify the full path to the library in the input file:: + + thermoLibraries = ['path/to/your/thermo/library/file.py'] + +or:: + + reactionLibraries = [(path/to/your/kinetic/library/folder/'] + +Combinations in any order of RMG's legacy libraries and users' external libraries are allowed, +and the order in which the libraries are specified is the order in which they are loaded and given priority. + .. _seedmechanism: Seed Mechanisms