From 34ed5a3d53350bdd88dbf1cd1a47c435b5863f65 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sat, 4 Sep 2021 18:36:23 -0400 Subject: [PATCH 01/24] add functionality for constant species in Reactor classes --- rmgpy/rmg/reactors.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/rmgpy/rmg/reactors.py b/rmgpy/rmg/reactors.py index 522cb33af4..d102d36c32 100644 --- a/rmgpy/rmg/reactors.py +++ b/rmgpy/rmg/reactors.py @@ -365,6 +365,10 @@ def __init__(self, core_phase_system, edge_phase_system, initial_conditions, ter self.terminations = [to_rms(term) for term in terminations] + def get_const_spc_indices(self,spcs=None): + rms_species_names = [spc.name for spc in self.core_phase_system.get_rms_species_list()] + return [rms_species_names.index(name) for name in self.const_spc_names] + def finish_termination_criteria(self): """ Convert tuples into TerminationConversion objects From a6742e9b7353f51708828da402dbb607fbf52dc5 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sat, 4 Sep 2021 18:41:28 -0400 Subject: [PATCH 02/24] For surface species calculate solvent correction from the saturated desorbed species --- rmgpy/data/solvation.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/rmgpy/data/solvation.py b/rmgpy/data/solvation.py index e961e503ef..99a5e26979 100644 --- a/rmgpy/data/solvation.py +++ b/rmgpy/data/solvation.py @@ -1198,6 +1198,9 @@ def get_solute_data_from_groups(self, species): by gas-phase thermo estimate. """ molecule = species.molecule[0] + if molecule.contains_surface_site(): + molecule = molecule.get_desorbed_molecules()[0] + molecule.saturate_unfilled_valence() molecule.clear_labeled_atoms() molecule.update_atomtypes() solute_data = self.estimate_solute_via_group_additivity(molecule) From 49bb3c102320117efa26acb678d1701f6fa14c28 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 2 Sep 2021 14:49:18 -0400 Subject: [PATCH 03/24] added `max_surface_bond_order` constraint --- rmgpy/constraints.py | 6 ++++++ rmgpy/rmg/input.py | 1 + 2 files changed, 7 insertions(+) diff --git a/rmgpy/constraints.py b/rmgpy/constraints.py index 14de4275bc..30011c381b 100644 --- a/rmgpy/constraints.py +++ b/rmgpy/constraints.py @@ -113,6 +113,12 @@ def fails_species_constraints(species): if struct.get_num_atoms('X') > max_surface_sites: return True + max_surface_bond_order = species_constraints.get('maximumSurfaceBondOrder', -1) + if max_surface_sites != -1: + for site in struct.get_surface_sites(): + if site.get_total_bond_order() > max_surface_bond_order: + return True + max_heavy_atoms = species_constraints.get('maximumHeavyAtoms', -1) if max_heavy_atoms != -1: if struct.get_num_atoms() - struct.get_num_atoms('H') > max_heavy_atoms: diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index 9138b887af..271513f503 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -1374,6 +1374,7 @@ def generated_species_constraints(**kwargs): 'maximumSiliconAtoms', 'maximumSulfurAtoms', 'maximumSurfaceSites', + 'maximumSurfaceBondOrder', 'maximumHeavyAtoms', 'maximumRadicalElectrons', 'maximumSingletCarbenes', From 7cee1425d1ca7e7f85783e45bc4c032907510f77 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Mon, 4 Oct 2021 17:44:28 -0400 Subject: [PATCH 04/24] fix bug with max_surface_sites specification --- rmgpy/constraints.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/constraints.py b/rmgpy/constraints.py index 30011c381b..8f9ac10638 100644 --- a/rmgpy/constraints.py +++ b/rmgpy/constraints.py @@ -114,7 +114,7 @@ def fails_species_constraints(species): return True max_surface_bond_order = species_constraints.get('maximumSurfaceBondOrder', -1) - if max_surface_sites != -1: + if max_surface_bond_order != -1: for site in struct.get_surface_sites(): if site.get_total_bond_order() > max_surface_bond_order: return True From 5badbd1950d72f81ae2e537b3a7beed5009e1fc6 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Mon, 4 Oct 2021 17:42:45 -0400 Subject: [PATCH 05/24] handle case when symbol is in all caps --- arkane/common.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/arkane/common.py b/arkane/common.py index 4623a821ac..3c50c47699 100644 --- a/arkane/common.py +++ b/arkane/common.py @@ -390,7 +390,11 @@ def get_element_mass(input_element, isotope=None): number = input_element elif isinstance(input_element, str): symbol = input_element - number = next(key for key, value in symbol_by_number.items() if value == input_element) + try: + number = next(key for key, value in symbol_by_number.items() if value == input_element) + except: + symbol = input_element[0] + input_element[1].lower() + number = [key for key, value in symbol_by_number.items() if value == symbol][0] if symbol is None or number is None: raise ValueError('Could not identify element {0}'.format(input_element)) From 92fad47652baa92029f97c56ec8f7294ed70560a Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Mon, 4 Oct 2021 17:48:25 -0400 Subject: [PATCH 06/24] make it possible to explicitly forbid molecules and groups in the input file --- rmgpy/constraints.py | 10 ++++++++++ rmgpy/rmg/input.py | 8 +++++--- rmgpy/rmg/main.py | 2 +- 3 files changed, 16 insertions(+), 4 deletions(-) diff --git a/rmgpy/constraints.py b/rmgpy/constraints.py index 8f9ac10638..83c1151795 100644 --- a/rmgpy/constraints.py +++ b/rmgpy/constraints.py @@ -83,6 +83,16 @@ def fails_species_constraints(species): if struct.is_isomorphic(molecule): return False + explicitly_forbidden_molecules = species_constraints.get('explicitlyForbiddenMolecules', []) + for molecule in explicitly_forbidden_molecules: + if struct.is_isomorphic(molecule): + return True + + explicitly_forbidden_groups = species_constraints.get('explicitlyForbiddenGroups', []) + for group in explicitly_forbidden_groups: + if struct.is_subgraph_isomorphic(group): + return True + max_carbon_atoms = species_constraints.get('maximumCarbonAtoms', -1) if max_carbon_atoms != -1: if struct.get_num_atoms('C') > max_carbon_atoms: diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index 271513f503..c261eaf49b 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -36,8 +36,7 @@ from rmgpy import settings from rmgpy.data.base import Entry from rmgpy.exceptions import DatabaseError, InputError -from rmgpy.molecule import Molecule -from rmgpy.molecule.group import Group +from rmgpy.molecule import Molecule, Group from rmgpy.quantity import Quantity, Energy, RateCoefficient, SurfaceConcentration from rmgpy.rmg.model import CoreEdgeReactionModel from rmgpy.rmg.settings import ModelSettings, SimulatorSettings @@ -238,7 +237,7 @@ def inchi(string): def adjacency_list(string): return Molecule().from_adjacency_list(string) -def adjacency_list_group(string): +def group_adjacency_list(string): return Group().from_adjacency_list(string) def react(tups): @@ -1381,6 +1380,8 @@ def generated_species_constraints(**kwargs): 'maximumCarbeneRadicals', 'allowSingletO2', 'speciesCuttingThreshold', + 'explicitlyForbiddenMolecules', + 'explicitlyForbiddenGroups', ] for key, value in kwargs.items(): @@ -1528,6 +1529,7 @@ def read_input_file(path, rmg0): 'InChI': inchi, 'adjacencyList': adjacency_list, 'adjacencyListGroup': adjacency_list_group, + 'groupAdjacencyList': group_adjacency_list, 'react': react, 'simpleReactor': simple_reactor, 'constantVIdealGasReactor' : constant_V_ideal_gas_reactor, diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index 4af6a43256..2f659dc263 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -65,7 +65,7 @@ from rmgpy.kinetics.diffusionLimited import diffusion_limiter from rmgpy.data.vaporLiquidMassTransfer import vapor_liquid_mass_transfer from rmgpy.kinetics import ThirdBody -from rmgpy.molecule import Molecule +from rmgpy.molecule import Molecule, Group from rmgpy.qm.main import QMDatabaseWriter from rmgpy.reaction import Reaction from rmgpy.rmg.listener import SimulationProfileWriter, SimulationProfilePlotter From 54d9e72486f1346bdc5075e80c10b4f4eb402055 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sun, 7 May 2023 12:19:24 -0700 Subject: [PATCH 07/24] fig adjacencyListGroup in input --- rmgpy/rmg/input.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index c261eaf49b..c19fd8986f 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -1528,7 +1528,7 @@ def read_input_file(path, rmg0): 'SMILES': smiles, 'InChI': inchi, 'adjacencyList': adjacency_list, - 'adjacencyListGroup': adjacency_list_group, + 'adjacencyListGroup': group_adjacency_list, 'groupAdjacencyList': group_adjacency_list, 'react': react, 'simpleReactor': simple_reactor, From bd1d25dd154542acf5b83392ca755e49e5d8fffb Mon Sep 17 00:00:00 2001 From: Richard West Date: Fri, 14 Jul 2023 12:01:54 -0400 Subject: [PATCH 08/24] Revert "make it possible to explicitly forbid molecules and groups in the input file" This reverts commit 92fad47652baa92029f97c56ec8f7294ed70560a. also "fig adjacencyListGroup in input" which was a fixup. The functionality has been available for 2 years since https://github.com/ReactionMechanismGenerator/RMG-Py/pull/2185 so adding a new syntax for doing it would be confusing and unnecessarry. --- rmgpy/constraints.py | 10 ---------- rmgpy/rmg/input.py | 10 ++++------ rmgpy/rmg/main.py | 2 +- 3 files changed, 5 insertions(+), 17 deletions(-) diff --git a/rmgpy/constraints.py b/rmgpy/constraints.py index 83c1151795..8f9ac10638 100644 --- a/rmgpy/constraints.py +++ b/rmgpy/constraints.py @@ -83,16 +83,6 @@ def fails_species_constraints(species): if struct.is_isomorphic(molecule): return False - explicitly_forbidden_molecules = species_constraints.get('explicitlyForbiddenMolecules', []) - for molecule in explicitly_forbidden_molecules: - if struct.is_isomorphic(molecule): - return True - - explicitly_forbidden_groups = species_constraints.get('explicitlyForbiddenGroups', []) - for group in explicitly_forbidden_groups: - if struct.is_subgraph_isomorphic(group): - return True - max_carbon_atoms = species_constraints.get('maximumCarbonAtoms', -1) if max_carbon_atoms != -1: if struct.get_num_atoms('C') > max_carbon_atoms: diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index c19fd8986f..271513f503 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -36,7 +36,8 @@ from rmgpy import settings from rmgpy.data.base import Entry from rmgpy.exceptions import DatabaseError, InputError -from rmgpy.molecule import Molecule, Group +from rmgpy.molecule import Molecule +from rmgpy.molecule.group import Group from rmgpy.quantity import Quantity, Energy, RateCoefficient, SurfaceConcentration from rmgpy.rmg.model import CoreEdgeReactionModel from rmgpy.rmg.settings import ModelSettings, SimulatorSettings @@ -237,7 +238,7 @@ def inchi(string): def adjacency_list(string): return Molecule().from_adjacency_list(string) -def group_adjacency_list(string): +def adjacency_list_group(string): return Group().from_adjacency_list(string) def react(tups): @@ -1380,8 +1381,6 @@ def generated_species_constraints(**kwargs): 'maximumCarbeneRadicals', 'allowSingletO2', 'speciesCuttingThreshold', - 'explicitlyForbiddenMolecules', - 'explicitlyForbiddenGroups', ] for key, value in kwargs.items(): @@ -1528,8 +1527,7 @@ def read_input_file(path, rmg0): 'SMILES': smiles, 'InChI': inchi, 'adjacencyList': adjacency_list, - 'adjacencyListGroup': group_adjacency_list, - 'groupAdjacencyList': group_adjacency_list, + 'adjacencyListGroup': adjacency_list_group, 'react': react, 'simpleReactor': simple_reactor, 'constantVIdealGasReactor' : constant_V_ideal_gas_reactor, diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index 2f659dc263..4af6a43256 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -65,7 +65,7 @@ from rmgpy.kinetics.diffusionLimited import diffusion_limiter from rmgpy.data.vaporLiquidMassTransfer import vapor_liquid_mass_transfer from rmgpy.kinetics import ThirdBody -from rmgpy.molecule import Molecule, Group +from rmgpy.molecule import Molecule from rmgpy.qm.main import QMDatabaseWriter from rmgpy.reaction import Reaction from rmgpy.rmg.listener import SimulationProfileWriter, SimulationProfilePlotter From 20adbd110282575e3a047646b520bdd49e970db4 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 11 Jan 2022 13:39:36 -0500 Subject: [PATCH 09/24] add comment attribute to species and reactions in rms yaml file --- rmgpy/yml.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/rmgpy/yml.py b/rmgpy/yml.py index c21e57d7d5..bd41c45d25 100644 --- a/rmgpy/yml.py +++ b/rmgpy/yml.py @@ -125,6 +125,7 @@ def obj_to_dict(obj, spcs, names=None, label="solvent"): result_dict["henrylawconstant"]["type"] = "TemperatureDependentHenryLawConstant" result_dict["henrylawconstant"]["Ts"] = obj.henry_law_constant_data.Ts result_dict["henrylawconstant"]["kHs"] = obj.henry_law_constant_data.kHs + result_dict["comment"] = obj.thermo.comment elif isinstance(obj, NASA): result_dict["polys"] = [obj_to_dict(k, spcs) for k in obj.polynomials] result_dict["type"] = "NASA" @@ -140,6 +141,7 @@ def obj_to_dict(obj, spcs, names=None, label="solvent"): result_dict["type"] = "ElementaryReaction" result_dict["radicalchange"] = sum([get_radicals(x) for x in obj.products]) - \ sum([get_radicals(x) for x in obj.reactants]) + result_dict["comment"] = obj.kinetics.comment elif isinstance(obj, Arrhenius): obj.change_t0(1.0) result_dict["type"] = "Arrhenius" From 4425ac56446667a901057b6168b4665fc811f760 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 11 Jan 2022 14:13:28 -0500 Subject: [PATCH 10/24] when rule making is done on one processor don't use pool --- rmgpy/data/kinetics/family.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index a3d4e4b630..08e34482f3 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -3621,9 +3621,12 @@ def make_bm_rules_from_template_rxn_map(self, template_rxn_map, nprocs=1, Tref=1 inds = inds.tolist() revinds = [inds.index(x) for x in np.arange(len(inputs))] - pool = mp.Pool(nprocs) + if nprocs > 1: + pool = mp.Pool(nprocs) + kinetics_list = np.array(pool.map(_make_rule, inputs[inds])) + else: + kinetics_list = np.array(list(map(_make_rule, inputs[inds]))) - kinetics_list = np.array(pool.map(_make_rule, inputs[inds])) kinetics_list = kinetics_list[revinds] # fix order for i, kinetics in enumerate(kinetics_list): From 7d527a3003aa7ccd5277d4cb403c05d62cb3b463 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 11 Jan 2022 14:31:18 -0500 Subject: [PATCH 11/24] add function for averaging kinetics --- rmgpy/data/kinetics/family.py | 77 +++++++++++++++++++++++++++++++++++ 1 file changed, 77 insertions(+) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 08e34482f3..e49d352fba 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -4673,3 +4673,80 @@ def _child_make_tree_nodes(family, child_conn, template_rxn_map, obj, T, nprocs, extension_iter_max=extension_iter_max, extension_iter_item_cap=extension_iter_item_cap) child_conn.send(list(family.groups.entries.values())) + +def average_kinetics(kinetics_list): + """ + Based on averaging log k. + Hence we average n, Ea, arithmetically, but we + average log A (geometric average) + """ + logA = 0.0 + n = 0.0 + Ea = 0.0 + alpha = 0.5 + electrons = None + if isinstance(kinetics_list[0], SurfaceChargeTransfer) or isinstance(kinetics_list[0], ArrheniusChargeTransfer): + if electrons is None: + electrons = kinetics_list[0].electrons.value_si + assert all(np.abs(k.V0.value_si) < 0.0001 for k in kinetics_list), [k.V0.value_si for k in kinetics_list] + assert all(np.abs(k.alpha.value_si - 0.5) < 0.001 for k in kinetics_list), [k.alpha for k in kinetics_list] + V0 = 0.0 + count = 0 + for kinetics in kinetics_list: + count += 1 + logA += np.log10(kinetics.A.value_si) + n += kinetics.n.value_si + Ea += kinetics.Ea.value_si + + logA /= count + n /= count + alpha /= count + Ea /= count + Aunits = kinetics_list[0].A.units + if Aunits == 'cm^3/(mol*s)' or Aunits == 'cm^3/(molecule*s)' or Aunits == 'm^3/(molecule*s)': + Aunits = 'm^3/(mol*s)' + elif Aunits == 'cm^6/(mol^2*s)' or Aunits == 'cm^6/(molecule^2*s)' or Aunits == 'm^6/(molecule^2*s)': + Aunits = 'm^6/(mol^2*s)' + elif Aunits == 's^-1' or Aunits == 'm^3/(mol*s)' or Aunits == 'm^6/(mol^2*s)': + # they were already in SI + pass + elif Aunits in ['m^2/(mol*s)', 'cm^2/(mol*s)', 'm^2/(molecule*s)', 'cm^2/(molecule*s)']: + # surface: bimolecular (Langmuir-Hinshelwood) + Aunits = 'm^2/(mol*s)' + elif Aunits in ['m^5/(mol^2*s)', 'cm^5/(mol^2*s)', 'm^5/(molecule^2*s)', 'cm^5/(molecule^2*s)']: + # surface: dissociative adsorption + Aunits = 'm^5/(mol^2*s)' + elif Aunits == '': + # surface: sticking coefficient + pass + else: + raise Exception('Invalid units {0} for averaging kinetics.'.format(Aunits)) + + if type(kinetics) not in [Arrhenius,SurfaceChargeTransfer,ArrheniusChargeTransfer]: + raise Exception('Invalid kinetics type {0!r} for {1!r}.'.format(type(kinetics), self)) + + if isinstance(kinetics, SurfaceChargeTransfer): + averaged_kinetics = SurfaceChargeTransfer( + A=(10 ** logA, Aunits), + n=n, + electrons=electrons, + alpha=alpha, + V0=(V0,'V'), + Ea=(Ea * 0.001, "kJ/mol"), + ) + elif isinstance(kinetics, ArrheniusChargeTransfer): + averaged_kinetics = ArrheniusChargeTransfer( + A=(10 ** logA, Aunits), + n=n, + electrons=electrons, + alpha=alpha, + V0=(V0,'V'), + Ea=(Ea * 0.001, "kJ/mol"), + ) + else: + averaged_kinetics = Arrhenius( + A=(10 ** logA, Aunits), + n=n, + Ea=(Ea * 0.001, "kJ/mol"), + ) + return averaged_kinetics From 23d88f246451123c78835e118e235c3bf2efacc2 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sat, 18 Jun 2022 12:46:08 -0700 Subject: [PATCH 12/24] check reaction exists before checking coverage_dependence --- rmgpy/rmg/model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index 12c6ecdf80..3194ceec80 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -1713,7 +1713,7 @@ def add_reaction_library_to_edge(self, reaction_library): reversible=rxn.reversible ) r, isNew = self.make_new_reaction(rxn) # updates self.new_species_list and self.newReactionlist - if getattr(r.kinetics, 'coverage_dependence', None): + if r and getattr(r.kinetics, 'coverage_dependence', None): self.process_coverage_dependence(r.kinetics) if not isNew: logging.info("This library reaction was not new: {0}".format(rxn)) From 0a6c456c41424d3c279e22f1ab790d7ca0cd5ea2 Mon Sep 17 00:00:00 2001 From: Richard West Date: Sun, 2 Jul 2023 17:17:10 +0100 Subject: [PATCH 13/24] Replace Asserts with ValueErrors --- rmgpy/data/kinetics/family.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index e49d352fba..10c9c9a2d0 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -4688,8 +4688,10 @@ def average_kinetics(kinetics_list): if isinstance(kinetics_list[0], SurfaceChargeTransfer) or isinstance(kinetics_list[0], ArrheniusChargeTransfer): if electrons is None: electrons = kinetics_list[0].electrons.value_si - assert all(np.abs(k.V0.value_si) < 0.0001 for k in kinetics_list), [k.V0.value_si for k in kinetics_list] - assert all(np.abs(k.alpha.value_si - 0.5) < 0.001 for k in kinetics_list), [k.alpha for k in kinetics_list] + if not all(np.abs(k.V0.value_si) < 0.0001 for k in kinetics_list): + raise ValueError(f"Trying to average charge transfer rates with non-zero V0 values: {[k.V0.value_si for k in kinetics_list]}") + if not all(np.abs(k.alpha.value_si - 0.5) < 0.001 for k in kinetics_list): + raise ValueError(f"Trying to average charge transfer rates with alpha values not equal to 0.5: {[k.alpha for k in kinetics_list]}") V0 = 0.0 count = 0 for kinetics in kinetics_list: From 21801750f0396821f5700f0a6a344adecf3c927f Mon Sep 17 00:00:00 2001 From: Richard West Date: Sun, 2 Jul 2023 17:23:09 +0100 Subject: [PATCH 14/24] Remove charge transfer types from average_kinetics (to be reverted) SurfaceChargeTransfer and ArrheniusChargeTransfer are not yet defined (on this branch). Revert this commin once they are. --- rmgpy/data/kinetics/family.py | 33 +++------------------------------ 1 file changed, 3 insertions(+), 30 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 10c9c9a2d0..65d8571228 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -4683,16 +4683,6 @@ def average_kinetics(kinetics_list): logA = 0.0 n = 0.0 Ea = 0.0 - alpha = 0.5 - electrons = None - if isinstance(kinetics_list[0], SurfaceChargeTransfer) or isinstance(kinetics_list[0], ArrheniusChargeTransfer): - if electrons is None: - electrons = kinetics_list[0].electrons.value_si - if not all(np.abs(k.V0.value_si) < 0.0001 for k in kinetics_list): - raise ValueError(f"Trying to average charge transfer rates with non-zero V0 values: {[k.V0.value_si for k in kinetics_list]}") - if not all(np.abs(k.alpha.value_si - 0.5) < 0.001 for k in kinetics_list): - raise ValueError(f"Trying to average charge transfer rates with alpha values not equal to 0.5: {[k.alpha for k in kinetics_list]}") - V0 = 0.0 count = 0 for kinetics in kinetics_list: count += 1 @@ -4702,7 +4692,6 @@ def average_kinetics(kinetics_list): logA /= count n /= count - alpha /= count Ea /= count Aunits = kinetics_list[0].A.units if Aunits == 'cm^3/(mol*s)' or Aunits == 'cm^3/(molecule*s)' or Aunits == 'm^3/(molecule*s)': @@ -4724,27 +4713,11 @@ def average_kinetics(kinetics_list): else: raise Exception('Invalid units {0} for averaging kinetics.'.format(Aunits)) - if type(kinetics) not in [Arrhenius,SurfaceChargeTransfer,ArrheniusChargeTransfer]: + if type(kinetics) not in [Arrhenius,]: raise Exception('Invalid kinetics type {0!r} for {1!r}.'.format(type(kinetics), self)) - if isinstance(kinetics, SurfaceChargeTransfer): - averaged_kinetics = SurfaceChargeTransfer( - A=(10 ** logA, Aunits), - n=n, - electrons=electrons, - alpha=alpha, - V0=(V0,'V'), - Ea=(Ea * 0.001, "kJ/mol"), - ) - elif isinstance(kinetics, ArrheniusChargeTransfer): - averaged_kinetics = ArrheniusChargeTransfer( - A=(10 ** logA, Aunits), - n=n, - electrons=electrons, - alpha=alpha, - V0=(V0,'V'), - Ea=(Ea * 0.001, "kJ/mol"), - ) + if False: + pass else: averaged_kinetics = Arrhenius( A=(10 ** logA, Aunits), From 9609275d46511b7bb8c8d5b32bd5ca86f551f37a Mon Sep 17 00:00:00 2001 From: Richard West Date: Tue, 4 Jul 2023 18:33:01 +0100 Subject: [PATCH 15/24] Unit test for average_kinetics method For now it just tests Arrhenius types, but should eventually work on others. --- rmgpy/data/kinetics/familyTest.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/rmgpy/data/kinetics/familyTest.py b/rmgpy/data/kinetics/familyTest.py index af39d28659..d3b361bd40 100644 --- a/rmgpy/data/kinetics/familyTest.py +++ b/rmgpy/data/kinetics/familyTest.py @@ -37,12 +37,15 @@ import numpy as np from rmgpy import settings +import rmgpy.data.kinetics.family from rmgpy.data.kinetics.database import KineticsDatabase from rmgpy.data.kinetics.family import TemplateReaction from rmgpy.data.rmg import RMGDatabase from rmgpy.data.thermo import ThermoDatabase from rmgpy.molecule import Molecule from rmgpy.species import Species +from rmgpy.kinetics import Arrhenius + ################################################### @@ -1090,6 +1093,23 @@ def test_retaining_atom_labels_in_template_reaction(self): self.assertEqual([(label, str(atom)) for label, atom in reaction_list_2[0].labeled_atoms['products'].items()], [('*1', 'C'), ('*2', 'C.'), ('*3', 'H')]) + + def test_average_kinetics(self): + """ + Test that the average kinetics are calculated correctly + """ + k1 = Arrhenius(A=(1e+13, 'cm^3/(mol*s)'), n=0, Ea=(0, 'kJ/mol'), T0=(1, 'K')) + k2 = Arrhenius(A=(4e+13, 'cm^3/(mol*s)'), n=1, Ea=(10, 'kJ/mol'), T0=(1, 'K')) + kav = rmgpy.data.kinetics.family.average_kinetics([k1, k2]) + self.assertAlmostEqual(kav.A.value_si, 2.0e7, 2) # m3/mol/s + self.assertAlmostEqual(kav.n.value_si, 0.5, 6) + self.assertAlmostEqual(kav.Ea.value_si, 5.0e3, 2) + self.assertAlmostEqual(kav.T0.value_si, 1.0, 6) + self.assertEqual(kav.A.units, 'm^3/(mol*s)') + self.assertAlmostEqual(np.log(kav.get_rate_coefficient(300)), + np.average([np.log(k1.get_rate_coefficient(300)), + np.log(k2.get_rate_coefficient(300))]), 6) + ################################################################################ From 88adf99e9e670377d85a6f33dd8172eb89e95a0c Mon Sep 17 00:00:00 2001 From: Richard West Date: Fri, 14 Jul 2023 11:22:24 -0400 Subject: [PATCH 16/24] Add documentation (examples) for recently added species constraints. maximumSurfaceSites and maximumSurfaceBondOrder --- documentation/source/users/rmg/input.rst | 1 + examples/rmg/commented/input.py | 2 ++ 2 files changed, 3 insertions(+) diff --git a/documentation/source/users/rmg/input.rst b/documentation/source/users/rmg/input.rst index 7b9989460d..33814e3f84 100644 --- a/documentation/source/users/rmg/input.rst +++ b/documentation/source/users/rmg/input.rst @@ -989,6 +989,7 @@ all of RMG's reaction families. :: maximumSulfurAtoms=2, maximumHeavyAtoms=10, maximumSurfaceSites=2, + maximumSurfaceBondOrder=4, maximumRadicalElectrons=2, maximumSingletCarbenes=1, maximumCarbeneRadicals=0, diff --git a/examples/rmg/commented/input.py b/examples/rmg/commented/input.py index 4ae3a23e6b..91b01b2b11 100644 --- a/examples/rmg/commented/input.py +++ b/examples/rmg/commented/input.py @@ -271,6 +271,8 @@ maximumNitrogenAtoms=0, maximumSiliconAtoms=0, maximumSulfurAtoms=0, + maximumSurfaceSites=2, # maximum number of surface sites (for heterogeneous catalysis) + maximumSurfaceBondOrder=2, # maximum bond order of each surface sites (for heterogeneous catalysis) # max number of non-hydrogen atoms # maximumHeavyAtoms=20, # maximum radicals on a molecule From f8fdf43d9f982990a42b046799382980f609b977 Mon Sep 17 00:00:00 2001 From: Richard West Date: Fri, 14 Jul 2023 11:23:24 -0400 Subject: [PATCH 17/24] Check species constraints in consistent order. Doesn't make a real difference, just easier for a human to parse the code without making errors, when the "if" blocks are in the same order as the documentation and examples. --- rmgpy/constraints.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/rmgpy/constraints.py b/rmgpy/constraints.py index 8f9ac10638..bdd6607f9f 100644 --- a/rmgpy/constraints.py +++ b/rmgpy/constraints.py @@ -108,6 +108,11 @@ def fails_species_constraints(species): if struct.get_num_atoms('S') > max_sulfur_atoms: return True + max_heavy_atoms = species_constraints.get('maximumHeavyAtoms', -1) + if max_heavy_atoms != -1: + if struct.get_num_atoms() - struct.get_num_atoms('H') > max_heavy_atoms: + return True + max_surface_sites = species_constraints.get('maximumSurfaceSites', -1) if max_surface_sites != -1: if struct.get_num_atoms('X') > max_surface_sites: @@ -119,11 +124,6 @@ def fails_species_constraints(species): if site.get_total_bond_order() > max_surface_bond_order: return True - max_heavy_atoms = species_constraints.get('maximumHeavyAtoms', -1) - if max_heavy_atoms != -1: - if struct.get_num_atoms() - struct.get_num_atoms('H') > max_heavy_atoms: - return True - max_radicals = species_constraints.get('maximumRadicalElectrons', -1) if max_radicals != -1: if struct.get_radical_count() > max_radicals: From f9991fee9626f793ae171d11c9eff43e3be43dfe Mon Sep 17 00:00:00 2001 From: Richard West Date: Fri, 14 Jul 2023 11:23:43 -0400 Subject: [PATCH 18/24] Unit test for maximumSurfaceBondOrder species constraint. --- rmgpy/constraintsTest.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/rmgpy/constraintsTest.py b/rmgpy/constraintsTest.py index 47bceeebdb..04e6770ad9 100644 --- a/rmgpy/constraintsTest.py +++ b/rmgpy/constraintsTest.py @@ -62,6 +62,7 @@ def setUpClass(cls): maximumSiliconAtoms=1, maximumSulfurAtoms=1, maximumSurfaceSites=2, + maximumSurfaceBondOrder=3, maximumHeavyAtoms=3, maximumRadicalElectrons=2, maximumSingletCarbenes=1, @@ -211,6 +212,16 @@ def test_surface_site_constraint(self): self.rmg.species_constraints['maximumCarbonAtoms'] = max_carbon self.rmg.species_constraints['maximumHeavyAtoms'] = max_heavy_atoms + def test_surface_bond_order_constraint(self): + """ + Test that we can constrain the max bond order of surface sites. + """ + mol_1site = Molecule().from_adjacency_list(""" +1 C u0 p0 c0 {2,Q} +2 X u0 p0 c0 {1,Q} +""") + self.assertTrue(fails_species_constraints(mol_1site)) + def test_heavy_constraint(self): """ Test that we can constrain the max number of heavy atoms. From 25ff54899b28458f16a6da4ae6d5f4466f517afe Mon Sep 17 00:00:00 2001 From: Richard West Date: Fri, 14 Jul 2023 14:17:33 -0400 Subject: [PATCH 19/24] Changing 'or or' to 'in set' checks. If you happen to pick the first of the list then the or or method was faster. But if you pick the last of the list (or not in the list) then the in tuple method is faster. Also changed some lists to sets for consistency, and because they're a touch faster. See pylint https://pylint.pycqa.org/en/latest/user_guide/messages/refactor/consider-using-in.html (Change requested in code review) Also switched a format to an fstring. --- rmgpy/data/kinetics/family.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 65d8571228..528b3f8da4 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -4694,27 +4694,27 @@ def average_kinetics(kinetics_list): n /= count Ea /= count Aunits = kinetics_list[0].A.units - if Aunits == 'cm^3/(mol*s)' or Aunits == 'cm^3/(molecule*s)' or Aunits == 'm^3/(molecule*s)': + if Aunits in {'cm^3/(mol*s)', 'cm^3/(molecule*s)', 'm^3/(molecule*s)'}: Aunits = 'm^3/(mol*s)' - elif Aunits == 'cm^6/(mol^2*s)' or Aunits == 'cm^6/(molecule^2*s)' or Aunits == 'm^6/(molecule^2*s)': + elif Aunits in {'cm^6/(mol^2*s)', 'cm^6/(molecule^2*s)', 'm^6/(molecule^2*s)'}: Aunits = 'm^6/(mol^2*s)' - elif Aunits == 's^-1' or Aunits == 'm^3/(mol*s)' or Aunits == 'm^6/(mol^2*s)': + elif Aunits in {'s^-1', 'm^3/(mol*s)', 'm^6/(mol^2*s)'}: # they were already in SI pass - elif Aunits in ['m^2/(mol*s)', 'cm^2/(mol*s)', 'm^2/(molecule*s)', 'cm^2/(molecule*s)']: + elif Aunits in {'m^2/(mol*s)', 'cm^2/(mol*s)', 'm^2/(molecule*s)', 'cm^2/(molecule*s)'}: # surface: bimolecular (Langmuir-Hinshelwood) Aunits = 'm^2/(mol*s)' - elif Aunits in ['m^5/(mol^2*s)', 'cm^5/(mol^2*s)', 'm^5/(molecule^2*s)', 'cm^5/(molecule^2*s)']: + elif Aunits in {'m^5/(mol^2*s)', 'cm^5/(mol^2*s)', 'm^5/(molecule^2*s)', 'cm^5/(molecule^2*s)'}: # surface: dissociative adsorption Aunits = 'm^5/(mol^2*s)' elif Aunits == '': # surface: sticking coefficient pass else: - raise Exception('Invalid units {0} for averaging kinetics.'.format(Aunits)) + raise Exception(f'Invalid units {Aunits} for averaging kinetics.') if type(kinetics) not in [Arrhenius,]: - raise Exception('Invalid kinetics type {0!r} for {1!r}.'.format(type(kinetics), self)) + raise Exception(f'Invalid kinetics type {type(kinetics)!r} for {self!r}.') if False: pass From 409b39dc5b24757746146fad19a7baaf3672bbc8 Mon Sep 17 00:00:00 2001 From: Richard West Date: Fri, 14 Jul 2023 16:48:08 -0400 Subject: [PATCH 20/24] Remove obsolete comments. They describe some kinetics comments that were moved in 921d425ac6083988290754381c3d6f75af8f420b and removed in 22e6a1352a830d1317664e02d4d9bea1ee590a2f --- rmgpy/rmg/model.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index 3194ceec80..e69c432ea8 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -1760,9 +1760,6 @@ def add_reaction_library_to_edge(self, reaction_library): self.add_species_to_edge(spec) for rxn in self.new_reaction_list: - # Note that we haven't actually evaluated any fluxes at this point - # Instead, we remove the comment below if the reaction is moved to - # the core later in the mechanism generation if not (self.pressure_dependence and rxn.elementary_high_p and rxn.is_unimolecular() and isinstance(rxn, LibraryReaction) and isinstance(rxn.kinetics, Arrhenius) and \ (self.pressure_dependence.maximum_atoms is None or self.pressure_dependence.maximum_atoms >= \ From 6071531d34c2d38c97ec0e5b7535e389b3357d95 Mon Sep 17 00:00:00 2001 From: Richard West Date: Fri, 14 Jul 2023 16:49:04 -0400 Subject: [PATCH 21/24] Update attribute name in a comment. Must be left over from the PEP8-ify project. --- rmgpy/rmg/model.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index e69c432ea8..b372421574 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -1608,7 +1608,7 @@ def add_seed_mechanism_to_core(self, seed_mechanism, react=False): kinetics=rxn.kinetics, duplicate=rxn.duplicate, reversible=rxn.reversible ) - r, isNew = self.make_new_reaction(rxn) # updates self.new_species_list and self.newReactionlist + r, isNew = self.make_new_reaction(rxn) # updates self.new_species_list and self.new_reaction_list if getattr(r.kinetics, 'coverage_dependence', None): self.process_coverage_dependence(r.kinetics) if not isNew: @@ -1712,8 +1712,8 @@ def add_reaction_library_to_edge(self, reaction_library): kinetics=rxn.kinetics, duplicate=rxn.duplicate, reversible=rxn.reversible ) - r, isNew = self.make_new_reaction(rxn) # updates self.new_species_list and self.newReactionlist if r and getattr(r.kinetics, 'coverage_dependence', None): + r, isNew = self.make_new_reaction(rxn) # updates self.new_species_list and self.new_reaction_list self.process_coverage_dependence(r.kinetics) if not isNew: logging.info("This library reaction was not new: {0}".format(rxn)) From 09a39e9e892fe8a1eb9b24dcd65295aab1270431 Mon Sep 17 00:00:00 2001 From: Richard West Date: Fri, 14 Jul 2023 16:51:09 -0400 Subject: [PATCH 22/24] Make check explicitly 'not None'. The only way r would be 'None' is if a symmetrical reaction (reactants == products) was found in a library, and hit this check in check_for_existing_reaction that's called by make_new_reaction. But that suggests some other problem. So I doubt this test is often needed. But Matt added it for something. I can't see any way it would be False, so now we check it's not None. (as requested in code review) --- rmgpy/rmg/model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index b372421574..b4134f60a0 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -1712,8 +1712,8 @@ def add_reaction_library_to_edge(self, reaction_library): kinetics=rxn.kinetics, duplicate=rxn.duplicate, reversible=rxn.reversible ) - if r and getattr(r.kinetics, 'coverage_dependence', None): r, isNew = self.make_new_reaction(rxn) # updates self.new_species_list and self.new_reaction_list + if r is not None and getattr(rxn.kinetics, 'coverage_dependence', None): self.process_coverage_dependence(r.kinetics) if not isNew: logging.info("This library reaction was not new: {0}".format(rxn)) From 4013707f24b54e788d2c339c5e5c68bc06efaabe Mon Sep 17 00:00:00 2001 From: Richard West Date: Fri, 14 Jul 2023 17:18:32 -0400 Subject: [PATCH 23/24] Comment an elegant alternative averaging scheme. Unfortunately we need to revert the commit titled "Remove charge transfer types from average_kinetics (to be reverted)" at some point in the not so distant future, and doing this different averaging scheme now would mess that up. --- rmgpy/data/kinetics/family.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 528b3f8da4..ff81d561c7 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -4693,6 +4693,12 @@ def average_kinetics(kinetics_list): logA /= count n /= count Ea /= count + + ## The above could be replaced with something like: + # logA, n, Ea = np.mean([[np.log10(k.A.value_si), + # k.n.value_si, + # k.Ea.value_si] for k in kinetics_list], axis=1) + Aunits = kinetics_list[0].A.units if Aunits in {'cm^3/(mol*s)', 'cm^3/(molecule*s)', 'm^3/(molecule*s)'}: Aunits = 'm^3/(mol*s)' From cfe970f4bcaf4037696c7625f1439811c3db678f Mon Sep 17 00:00:00 2001 From: Richard West Date: Fri, 14 Jul 2023 18:09:33 -0400 Subject: [PATCH 24/24] Simplify the atomic symbol to number conversion. --- arkane/common.py | 9 +++++---- arkane/commonTest.py | 1 + 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/arkane/common.py b/arkane/common.py index 3c50c47699..4f67ff1977 100644 --- a/arkane/common.py +++ b/arkane/common.py @@ -391,10 +391,10 @@ def get_element_mass(input_element, isotope=None): elif isinstance(input_element, str): symbol = input_element try: - number = next(key for key, value in symbol_by_number.items() if value == input_element) - except: - symbol = input_element[0] + input_element[1].lower() - number = [key for key, value in symbol_by_number.items() if value == symbol][0] + number = number_by_symbol[symbol] + except KeyError: + symbol = input_element.capitalize() + number = number_by_symbol[symbol] if symbol is None or number is None: raise ValueError('Could not identify element {0}'.format(input_element)) @@ -438,6 +438,7 @@ def get_element_mass(input_element, isotope=None): 92: 'U', 93: 'Np', 94: 'Pu', 95: 'Am', 96: 'Cm', 97: 'Bk', 98: 'Cf', 99: 'Es', 100: 'Fm', 101: 'Md', 102: 'No', 103: 'Lr', 104: 'Rf', 105: 'Db', 106: 'Sg', 107: 'Bh', 108: 'Hs', 109: 'Mt', 110: 'Ds', 111: 'Rg', 112: 'Cn', 113: 'Nh', 114: 'Fl', 115: 'Mc', 116: 'Lv', 117: 'Ts', 118: 'Og'} +number_by_symbol = {value: key for key, value in symbol_by_number.items()} # Structure of mass_by_symbol items: list(list(isotope1, mass1, weight1), list(isotope2, mass2, weight2), ...) mass_by_symbol = { diff --git a/arkane/commonTest.py b/arkane/commonTest.py index 7bca459d7a..50bff7fe77 100644 --- a/arkane/commonTest.py +++ b/arkane/commonTest.py @@ -459,6 +459,7 @@ def test_get_mass(self): """Test that the correct mass/number/isotope is returned from get_element_mass""" self.assertEqual(get_element_mass(1), (1.00782503224, 1)) # test input by integer self.assertEqual(get_element_mass('Si'), (27.97692653465, 14)) # test string input and most common isotope + self.assertEqual(get_element_mass('SI'), (27.97692653465, 14)) # test string in all caps self.assertEqual(get_element_mass('C', 13), (13.00335483507, 6)) # test specific isotope self.assertEqual(get_element_mass('Bk'), (247.0703073, 97)) # test a two-element array (no isotope data)