From 1065c9808a05c7f6051c607a5bca8353ec87d01e Mon Sep 17 00:00:00 2001 From: jonwzheng Date: Tue, 8 Jul 2025 15:57:52 -0400 Subject: [PATCH 01/10] update gri_mech_rxn_lib with updated constraints --- examples/rmg/gri_mech_rxn_lib/input.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/examples/rmg/gri_mech_rxn_lib/input.py b/examples/rmg/gri_mech_rxn_lib/input.py index a4749408099..033663936d6 100644 --- a/examples/rmg/gri_mech_rxn_lib/input.py +++ b/examples/rmg/gri_mech_rxn_lib/input.py @@ -40,6 +40,10 @@ maximumEdgeSpecies=100000 ) +generatedSpeciesConstraints( + maximumCarbeneRadicals=2, +) + options( units='si', generateOutputHTML=False, From 6794afd7563a0c2cb153381cb0ea1edfdb68d83e Mon Sep 17 00:00:00 2001 From: jonwzheng Date: Tue, 8 Jul 2025 16:19:50 -0400 Subject: [PATCH 02/10] Make constraints more descriptive Previously, failing a species constraint would not give the user any specific reason why the species failed the constraint check. This refactors the fails_species_constraint code so that it also returns a reason for the test failure. This required a small refactor to existing usages of the function. --- rmgpy/constraints.py | 32 ++++++------- rmgpy/data/kinetics/family.py | 9 +++- rmgpy/rmg/main.py | 5 ++- rmgpy/rmg/model.py | 10 +++-- test/rmgpy/constraintsTest.py | 84 +++++++++++++++++++++++------------ 5 files changed, 89 insertions(+), 51 deletions(-) diff --git a/rmgpy/constraints.py b/rmgpy/constraints.py index d00d1ca43bf..8a729769626 100644 --- a/rmgpy/constraints.py +++ b/rmgpy/constraints.py @@ -61,7 +61,8 @@ def pass_cutting_threshold(species): def fails_species_constraints(species): """ Pass in either a `Species` or `Molecule` object and checks whether it passes - the speciesConstraints set by the user. If not, returns `True` for failing speciesConstraints. + the speciesConstraints set by the user. If not, returns a tuple `(True, reason)` for failing speciesConstraints, + where `reason` is a string describing which constraint failed. If all constraints pass, returns `(False, None)`. """ from rmgpy.rmg.input import get_input @@ -81,62 +82,63 @@ def fails_species_constraints(species): explicitly_allowed_molecules = species_constraints.get('explicitlyAllowedMolecules', []) for molecule in explicitly_allowed_molecules: if struct.is_isomorphic(molecule): - return False + return (False, None) max_carbon_atoms = species_constraints.get('maximumCarbonAtoms', -1) if max_carbon_atoms != -1: if struct.get_num_atoms('C') > max_carbon_atoms: - return True + return (True, f"Exceeded maximumCarbonAtoms: {struct.get_num_atoms('C')} > {max_carbon_atoms}") max_oxygen_atoms = species_constraints.get('maximumOxygenAtoms', -1) if max_oxygen_atoms != -1: if struct.get_num_atoms('O') > max_oxygen_atoms: - return True + return (True, f"Exceeded maximumOxygenAtoms: {struct.get_num_atoms('O')} > {max_oxygen_atoms}") max_nitrogen_atoms = species_constraints.get('maximumNitrogenAtoms', -1) if max_nitrogen_atoms != -1: if struct.get_num_atoms('N') > max_nitrogen_atoms: - return True + return (True, f"Exceeded maximumNitrogenAtoms: {struct.get_num_atoms('N')} > {max_nitrogen_atoms}") max_silicon_atoms = species_constraints.get('maximumSiliconAtoms', -1) if max_silicon_atoms != -1: if struct.get_num_atoms('Si') > max_silicon_atoms: - return True + return (True, f"Exceeded maximumSiliconAtoms: {struct.get_num_atoms('Si')} > {max_silicon_atoms}") max_sulfur_atoms = species_constraints.get('maximumSulfurAtoms', -1) if max_sulfur_atoms != -1: if struct.get_num_atoms('S') > max_sulfur_atoms: - return True + return (True, f"Exceeded maximumSulfurAtoms: {struct.get_num_atoms('S')} > {max_sulfur_atoms}") 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 + heavy_atoms = struct.get_num_atoms() - struct.get_num_atoms('H') + if heavy_atoms > max_heavy_atoms: + return (True, f"Exceeded maximumHeavyAtoms: {heavy_atoms} > {max_heavy_atoms}") max_surface_sites = species_constraints.get('maximumSurfaceSites', -1) if max_surface_sites != -1: if struct.get_num_atoms('X') > max_surface_sites: - return True + return (True, f"Exceeded maximumSurfaceSites: {struct.get_num_atoms('X')} > {max_surface_sites}") max_surface_bond_order = species_constraints.get('maximumSurfaceBondOrder', -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 + return (True, f"Exceeded maximumSurfaceBondOrder at site: {site.get_total_bond_order()} > {max_surface_bond_order}") max_radicals = species_constraints.get('maximumRadicalElectrons', -1) if max_radicals != -1: if struct.get_radical_count() > max_radicals: - return True + return (True, f"Exceeded maximumRadicalElectrons: {struct.get_radical_count()} > {max_radicals}") max_carbenes = species_constraints.get('maximumSingletCarbenes', 1) if max_radicals != -1: if struct.get_singlet_carbene_count() > max_carbenes: - return True + return (True, f"Exceeded maximumSingletCarbenes: {struct.get_singlet_carbene_count()} > {max_carbenes}") max_carbene_radicals = species_constraints.get('maximumCarbeneRadicals', 0) if max_carbene_radicals != -1: if struct.get_singlet_carbene_count() > 0 and struct.get_radical_count() > max_carbene_radicals: - return True + return (True, f"Exceeded maximumCarbeneRadicals: {struct.get_radical_count()} > {max_carbene_radicals}") - return False + return (False, None) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 4bb0ee949e7..226bf646580 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -1654,8 +1654,13 @@ def _generate_product_structures(self, reactant_structures, maps, forward, relab for struct in product_structures: if self.is_molecule_forbidden(struct): raise ForbiddenStructureException() - if fails_species_constraints(struct): - raise ForbiddenStructureException() + failed, reason = fails_species_constraints(struct) + if failed: + raise ForbiddenStructureException( + "Species constraints forbids product species {0}. Please " + "reformulate constraints, or explicitly " + "allow it. Reason: {1}".format(struct.label, reason) + ) return product_structures diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index 88ea75ff4cc..e2b80967be1 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -689,14 +689,15 @@ def initialize(self, **kwargs): "Input species {0} is globally forbidden. You may explicitly " "allow it by adding 'input species' to the `generatedSpeciesConstraints` `allowed` list.".format(spec.label) ) - if fails_species_constraints(spec): + failed, reason = fails_species_constraints(spec) + if failed: if "allowed" in self.species_constraints and "input species" in self.species_constraints["allowed"]: self.species_constraints["explicitlyAllowedMolecules"].append(spec.molecule[0]) else: raise ForbiddenStructureException( "Species constraints forbids input species {0}. Please " "reformulate constraints, remove the species, or explicitly " - "allow it.".format(spec.label) + "allow it. Reason: {1}".format(spec.label, reason) ) # For liquidReactor, checks whether the solvent is listed as one of the initial species. diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index d677755b079..feaf92b5eb0 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -1673,14 +1673,15 @@ def add_seed_mechanism_to_core(self, seed_mechanism, react=False, requires_rms=F "found in a seed mechanism or reaction " "library.".format(spec.label, seed_mechanism.label) ) - if fails_species_constraints(spec): + failed, reason = fails_species_constraints(spec) + if failed: if "allowed" in rmg.species_constraints and "seed mechanisms" in rmg.species_constraints["allowed"]: rmg.species_constraints["explicitlyAllowedMolecules"].extend(spec.molecule) else: raise ForbiddenStructureException( "Species constraints forbids species {0} from seed mechanism {1}." " Please reformulate constraints, remove the species, or" - " explicitly allow it.".format(spec.label, seed_mechanism.label) + " explicitly allow it. Reason: {2}".format(spec.label, seed_mechanism.label, reason) ) for spec in edge_species_to_move+self.new_species_list: @@ -1800,14 +1801,15 @@ def add_reaction_library_to_edge(self, reaction_library, requires_rms=False): "inert unless found in a seed mechanism or reaction " "library.".format(spec.label, reaction_library.label) ) - if fails_species_constraints(spec): + failed, reason = fails_species_constraints(spec) + if failed: if "allowed" in rmg.species_constraints and "reaction libraries" in rmg.species_constraints["allowed"]: rmg.species_constraints["explicitlyAllowedMolecules"].extend(spec.molecule) else: raise ForbiddenStructureException( "Species constraints forbids species {0} from reaction library " "{1}. Please reformulate constraints, remove the species, or " - "explicitly allow it.".format(spec.label, reaction_library.label) + "explicitly allow it. Reason: {2}".format(spec.label, reaction_library.label, reason) ) for spec in self.new_species_list: diff --git a/test/rmgpy/constraintsTest.py b/test/rmgpy/constraintsTest.py index 45c411942c5..9db23f33a0c 100644 --- a/test/rmgpy/constraintsTest.py +++ b/test/rmgpy/constraintsTest.py @@ -83,8 +83,8 @@ def test_constraints_not_loaded(self, mock_logging): rmgpy.rmg.input.rmg = None mol = Molecule(smiles="C") - - assert not fails_species_constraints(mol) + failed, _ = fails_species_constraints(mol) + assert not failed mock_logging.debug.assert_called_with("Species constraints could not be found.") @@ -97,67 +97,80 @@ def test_species_input(self): """ spc = Species().from_smiles("C") - assert not fails_species_constraints(spc) + failed, _ = fails_species_constraints(spc) + assert not failed def test_explicitly_allowed_molecules(self): """ Test that we can explicitly allow molecules in species constraints. """ mol = Molecule(smiles="CCCC") - assert fails_species_constraints(mol) + failed, _ = fails_species_constraints(mol) + assert failed self.rmg.species_constraints["explicitlyAllowedMolecules"] = [Molecule(smiles="CCCC")] - assert not fails_species_constraints(mol) + failed, _ = fails_species_constraints(mol) + assert not failed def test_carbon_constraint(self): """ Test that we can constrain the max number of carbon atoms. """ mol1 = Molecule(smiles="CC") - assert not fails_species_constraints(mol1) + failed, _ = fails_species_constraints(mol1) + assert not failed mol2 = Molecule(smiles="CCC") - assert fails_species_constraints(mol2) + failed, _ = fails_species_constraints(mol2) + assert failed def test_oxygen_constraint(self): """ Test that we can constrain the max number of oxygen atoms. """ mol1 = Molecule(smiles="C=O") - assert not fails_species_constraints(mol1) + failed, _ = fails_species_constraints(mol1) + assert not failed mol2 = Molecule(smiles="OC=O") - assert fails_species_constraints(mol2) + failed, _ = fails_species_constraints(mol2) + assert failed def test_nitrogen_constraint(self): """ Test that we can constrain the max number of nitrogen atoms. """ mol1 = Molecule(smiles="CN") - assert not fails_species_constraints(mol1) + failed, _ = fails_species_constraints(mol1) + assert not failed mol2 = Molecule(smiles="NCN") - assert fails_species_constraints(mol2) + failed, _ = fails_species_constraints(mol2) + assert failed def test_silicon_constraint(self): """ Test that we can constrain the max number of silicon atoms. """ mol1 = Molecule(smiles="[SiH4]") - assert not fails_species_constraints(mol1) + failed, _ = fails_species_constraints(mol1) + assert not failed mol2 = Molecule(smiles="[SiH3][SiH3]") - assert fails_species_constraints(mol2) + failed, _ = fails_species_constraints(mol2) + assert failed def test_sulfur_constraint(self): """ Test that we can constrain the max number of sulfur atoms. """ mol1 = Molecule(smiles="CS") - assert not fails_species_constraints(mol1) + failed, _ = fails_species_constraints(mol1) + assert not failed mol2 = Molecule(smiles="SCS") - assert fails_species_constraints(mol2) + failed, _ = fails_species_constraints(mol2) + assert failed def test_surface_site_constraint(self): """ @@ -209,11 +222,17 @@ def test_surface_site_constraint(self): self.rmg.species_constraints["maximumCarbonAtoms"] = 3 self.rmg.species_constraints["maximumHeavyAtoms"] = 6 - assert not fails_species_constraints(mol_1site) - assert not fails_species_constraints(mol_2site) + failed, _ = fails_species_constraints(mol_1site) + assert not failed + + failed, _ = fails_species_constraints(mol_2site) + assert not failed + + failed, _ = fails_species_constraints(mol_3site_vdW) + assert failed - assert fails_species_constraints(mol_3site_vdW) - assert fails_species_constraints(mol_3site) + failed, _ = fails_species_constraints(mol_3site) + assert failed self.rmg.species_constraints["maximumCarbonAtoms"] = max_carbon self.rmg.species_constraints["maximumHeavyAtoms"] = max_heavy_atoms @@ -228,27 +247,32 @@ def test_surface_bond_order_constraint(self): 2 X u0 p0 c0 {1,Q} """ ) - assert fails_species_constraints(mol_1site) + failed, _ = fails_species_constraints(mol_1site) + assert failed def test_heavy_constraint(self): """ Test that we can constrain the max number of heavy atoms. """ mol1 = Molecule(smiles="CCO") - assert not fails_species_constraints(mol1) + failed, _ = fails_species_constraints(mol1) + assert not failed mol2 = Molecule(smiles="CCN=O") - assert fails_species_constraints(mol2) + failed, _ = fails_species_constraints(mol2) + assert failed def test_radical_constraint(self): """ Test that we can constrain the max number of radical electrons. """ mol1 = Molecule(smiles="[CH2][CH2]") - assert not fails_species_constraints(mol1) + failed, _ = fails_species_constraints(mol1) + assert not failed mol2 = Molecule(smiles="[CH2][CH][CH2]") - assert fails_species_constraints(mol2) + failed, _ = fails_species_constraints(mol2) + assert failed def test_carbene_constraint(self): """ @@ -261,7 +285,8 @@ def test_carbene_constraint(self): 3 H u0 p0 c0 {1,S} """ ) - assert not fails_species_constraints(mol1) + failed, _ = fails_species_constraints(mol1) + assert not failed mol2 = Molecule().from_adjacency_list( """ @@ -271,7 +296,8 @@ def test_carbene_constraint(self): 4 H u0 p0 c0 {3,S} """ ) - assert fails_species_constraints(mol2) + failed, _ = fails_species_constraints(mol2) + assert failed def test_carbene_radical_constraint(self): """ @@ -284,7 +310,8 @@ def test_carbene_radical_constraint(self): 3 H u0 p0 c0 {1,S} """ ) - assert not fails_species_constraints(mol1) + failed, _ = fails_species_constraints(mol1) + assert not failed mol2 = Molecule().from_adjacency_list( """ @@ -295,4 +322,5 @@ def test_carbene_radical_constraint(self): 5 H u0 p0 c0 {3,S} """ ) - assert fails_species_constraints(mol2) + failed, _ = fails_species_constraints(mol2) + assert failed From a41ed9a3889ec285e03f20bedec82d435f94115d Mon Sep 17 00:00:00 2001 From: jonwzheng Date: Tue, 8 Jul 2025 16:32:30 -0400 Subject: [PATCH 03/10] fix bug/typo with accessing nonexistent attribute --- rmgpy/data/kinetics/family.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 226bf646580..9f47187d037 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -1659,7 +1659,7 @@ def _generate_product_structures(self, reactant_structures, maps, forward, relab raise ForbiddenStructureException( "Species constraints forbids product species {0}. Please " "reformulate constraints, or explicitly " - "allow it. Reason: {1}".format(struct.label, reason) + "allow it. Reason: {1}".format(struct, reason) ) return product_structures From 9ee8bcc239dbe42ccdb95b69201ae67898062ec0 Mon Sep 17 00:00:00 2001 From: jonwzheng Date: Wed, 9 Jul 2025 12:38:14 -0400 Subject: [PATCH 04/10] update catalysis example path --- examples/rmg/catalysis/methane_steam/input.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/rmg/catalysis/methane_steam/input.py b/examples/rmg/catalysis/methane_steam/input.py index 22c4e098502..5eb8078b7d4 100644 --- a/examples/rmg/catalysis/methane_steam/input.py +++ b/examples/rmg/catalysis/methane_steam/input.py @@ -1,7 +1,7 @@ # Data sources database( thermoLibraries=['surfaceThermoNi111', 'surfaceThermoPt111', 'primaryThermoLibrary', 'thermo_DFT_CCSDTF12_BAC'], - reactionLibraries = [('Surface/Deutschmann_Ni', True)], # when Pt is used change the library to Surface/CPOX_Pt/Deutschmann2006 + reactionLibraries = [('Surface/Methane/Deutschmann_Ni', True)], # when Pt is used change the library to Surface/CPOX_Pt/Deutschmann2006 seedMechanisms = [], kineticsDepositories = ['training'], kineticsFamilies = ['surface','default'], From fa2aab9312b31915538cd9a964c084e2b7d2e316 Mon Sep 17 00:00:00 2001 From: jonwzheng Date: Wed, 9 Jul 2025 13:04:50 -0400 Subject: [PATCH 05/10] update seed mechanism loading kinetic libraries and families, which seem to be unhandled --- rmgpy/rmg/main.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index e2b80967be1..b3630846bac 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -587,9 +587,9 @@ def initialize(self, **kwargs): shutil.copyfile(self.species_map_path, os.path.join(filters_restart, "species_map.yml")) # Load the seed mechanism to get the core and edge species - self.database.kinetics.load_libraries(restart_dir, libraries=["restart", "restart_edge"]) + self.database.kinetics.load_libraries(restart_dir)#, libraries=["restart", "restart_edge"]) self.seed_mechanisms.append("restart") - self.reaction_libraries.append(("restart_edge", False)) +# self.reaction_libraries.append(("restart_edge", False)) # Set trimolecular reactant flags of reaction systems if self.trimolecular: From d82a7c5fe4b9ee921fdccb895647fe4432af048e Mon Sep 17 00:00:00 2001 From: jonwzheng Date: Wed, 9 Jul 2025 14:32:10 -0400 Subject: [PATCH 06/10] fix examples that don't exclude surface families --- examples/rmg/TEOS/input.py | 2 +- examples/rmg/e85/input.py | 2 +- examples/rmg/methylformate/input.py | 2 +- examples/rmg/minimal_sensitivity/input.py | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/rmg/TEOS/input.py b/examples/rmg/TEOS/input.py index 75b6a114085..86b3dcb77e8 100644 --- a/examples/rmg/TEOS/input.py +++ b/examples/rmg/TEOS/input.py @@ -4,7 +4,7 @@ reactionLibraries = [], seedMechanisms = [], kineticsDepositories = ['training'], - kineticsFamilies = ['!Intra_Disproportionation','!Substitution_O'], + kineticsFamilies = ['!Intra_Disproportionation','!Substitution_O', '!surface'], kineticsEstimator = 'rate rules', ) diff --git a/examples/rmg/e85/input.py b/examples/rmg/e85/input.py index 6dcf96aa433..66529b9280c 100644 --- a/examples/rmg/e85/input.py +++ b/examples/rmg/e85/input.py @@ -4,7 +4,7 @@ reactionLibraries = [], seedMechanisms = ['GRI-Mech3.0'], kineticsDepositories = ['training'], - kineticsFamilies = ['!Intra_Disproportionation','!Substitution_O'], + kineticsFamilies = ['!Intra_Disproportionation','!Substitution_O', '!surface'], kineticsEstimator = 'rate rules', ) diff --git a/examples/rmg/methylformate/input.py b/examples/rmg/methylformate/input.py index 783e04130e8..6cc01c85b1c 100644 --- a/examples/rmg/methylformate/input.py +++ b/examples/rmg/methylformate/input.py @@ -4,7 +4,7 @@ reactionLibraries = [('Methylformate',False),('Glarborg/highP',False)], seedMechanisms = ['Glarborg/C2'], kineticsDepositories = ['training'], - kineticsFamilies = ['!Intra_Disproportionation','!Substitution_O'], + kineticsFamilies = ['!Intra_Disproportionation','!Substitution_O', '!surface'], kineticsEstimator = 'rate rules', ) diff --git a/examples/rmg/minimal_sensitivity/input.py b/examples/rmg/minimal_sensitivity/input.py index c9358cadcd7..85eb1cdf73e 100644 --- a/examples/rmg/minimal_sensitivity/input.py +++ b/examples/rmg/minimal_sensitivity/input.py @@ -4,7 +4,7 @@ reactionLibraries = [], seedMechanisms = [], kineticsDepositories = ['training'], - kineticsFamilies = ['!Intra_Disproportionation','!Substitution_O'], + kineticsFamilies = ['!Intra_Disproportionation','!Substitution_O', '!surface'], kineticsEstimator = 'rate rules', ) From 2eb3abe96acb366f69cc6d747a9d0f8928927d2c Mon Sep 17 00:00:00 2001 From: jonwzheng Date: Wed, 9 Jul 2025 14:38:41 -0400 Subject: [PATCH 07/10] fix typo in raised error --- rmgpy/kinetics/chebyshev.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/kinetics/chebyshev.pyx b/rmgpy/kinetics/chebyshev.pyx index 3d0cb6bf1eb..7e8a24835f5 100644 --- a/rmgpy/kinetics/chebyshev.pyx +++ b/rmgpy/kinetics/chebyshev.pyx @@ -193,7 +193,7 @@ cdef class Chebyshev(PDepKineticsModel): raise KineticsError("The master equation data needs more temperature and pressure data " "points than are placed into Chebyshev polynomial. Currently, the " "data has {0} temperatures and the polynomial is set to have {1}. " - "The data has {2} pressures and the polynomial is set ot have {3}" + "The data has {2} pressures and the polynomial is set to have {3}" "".format(nT, degreeT, nP, degreeP)) elif nT < 1.25 * degreeT or nP < 1.25 * degreeP: logging.warning('This Chebyshev fitting has few degrees of freedom and may not be ' From 75631487ecea7bc9f7d0af705b8f3162dd612b51 Mon Sep 17 00:00:00 2001 From: jonwzheng Date: Wed, 9 Jul 2025 14:51:55 -0400 Subject: [PATCH 08/10] fix examples with malformed Chebyshev polynomials --- examples/rmg/MR_test/input.py | 4 ++-- examples/rmg/commented/input.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/rmg/MR_test/input.py b/examples/rmg/MR_test/input.py index 39f9e4d0fcf..6cbb04a2eb3 100644 --- a/examples/rmg/MR_test/input.py +++ b/examples/rmg/MR_test/input.py @@ -275,8 +275,8 @@ minimumNumberOfGrains=250, #the conditions for the rate to be output over #parameter order is: low_value, high_value, units, internal points - temperatures=(300,2200,'K',2), - pressures=(0.01,100.01,'bar',3), + temperatures=(300,2200,'K',10), + pressures=(0.01,100.01,'bar',10), #The two options for interpolation are 'PDepArrhenius' (no extra arguments) and #'Chebyshev' which is followed by the number of basis sets in #Temperature and Pressure. These values must be less than the number of diff --git a/examples/rmg/commented/input.py b/examples/rmg/commented/input.py index 0d6c78b5fea..e0044906c64 100644 --- a/examples/rmg/commented/input.py +++ b/examples/rmg/commented/input.py @@ -250,8 +250,8 @@ minimumNumberOfGrains=250, # the conditions for the rate to be output over # parameter order is: low_value, high_value, units, internal points - temperatures=(300, 2200, 'K', 2), - pressures=(0.01, 100, 'bar', 3), + temperatures=(300, 2200, 'K', 10), + pressures=(0.01, 100, 'bar', 10), # The two options for interpolation are 'PDepArrhenius' (no extra arguments) and # 'Chebyshev' which is followed by the number of basis sets in # Temperature and Pressure. These values must be less than the number of From 2f274dabf5602dd5d1927a811f671c25068c5597 Mon Sep 17 00:00:00 2001 From: jonwzheng Date: Wed, 9 Jul 2025 15:05:49 -0400 Subject: [PATCH 09/10] fix nox example excluding singlet O2 --- examples/rmg/nox_transitory_edge/input.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/rmg/nox_transitory_edge/input.py b/examples/rmg/nox_transitory_edge/input.py index a130c79fd06..569f116699c 100644 --- a/examples/rmg/nox_transitory_edge/input.py +++ b/examples/rmg/nox_transitory_edge/input.py @@ -12,10 +12,12 @@ reactive=True, structure=SMILES("CC"), ) + species( label='O2', structure=SMILES("[O][O]"), ) + species( label='N2', reactive=False, @@ -76,6 +78,7 @@ generatedSpeciesConstraints( allowed=['input species','seed mechanisms','reaction libraries'], + allowSingletO2=True, #maximumCarbonAtoms=5, #maximumOxygenAtoms=8, #maximumNitrogenAtoms=0, From 68bb109fc7e8b30edc5c3e4ad4d1fcebffc97368 Mon Sep 17 00:00:00 2001 From: jonwzheng Date: Fri, 2 Jan 2026 11:22:50 -0500 Subject: [PATCH 10/10] Update fails_species_constraints - returns only either False or a str, instead of a tuple - user must check if fails_species_constraints != False (or lazily if just fails_species_constraints) as both will return the reason for failing constraints; if no failure it will return False - Fixes a bug with implementation --- rmgpy/constraints.py | 32 +++++++------- rmgpy/data/kinetics/family.py | 4 +- rmgpy/rmg/main.py | 4 +- rmgpy/rmg/model.py | 8 ++-- test/rmgpy/constraintsTest.py | 81 ++++++++++++----------------------- 5 files changed, 51 insertions(+), 78 deletions(-) diff --git a/rmgpy/constraints.py b/rmgpy/constraints.py index 8a729769626..e88af0e3c02 100644 --- a/rmgpy/constraints.py +++ b/rmgpy/constraints.py @@ -61,8 +61,8 @@ def pass_cutting_threshold(species): def fails_species_constraints(species): """ Pass in either a `Species` or `Molecule` object and checks whether it passes - the speciesConstraints set by the user. If not, returns a tuple `(True, reason)` for failing speciesConstraints, - where `reason` is a string describing which constraint failed. If all constraints pass, returns `(False, None)`. + the speciesConstraints set by the user. If the species fails constraints, returns + a string `reason` describing which constraint failed. If all constraints pass, returns `False`. """ from rmgpy.rmg.input import get_input @@ -82,63 +82,63 @@ def fails_species_constraints(species): explicitly_allowed_molecules = species_constraints.get('explicitlyAllowedMolecules', []) for molecule in explicitly_allowed_molecules: if struct.is_isomorphic(molecule): - return (False, None) + return False max_carbon_atoms = species_constraints.get('maximumCarbonAtoms', -1) if max_carbon_atoms != -1: if struct.get_num_atoms('C') > max_carbon_atoms: - return (True, f"Exceeded maximumCarbonAtoms: {struct.get_num_atoms('C')} > {max_carbon_atoms}") + return f"Exceeded maximumCarbonAtoms: {struct.get_num_atoms('C')} > {max_carbon_atoms}" max_oxygen_atoms = species_constraints.get('maximumOxygenAtoms', -1) if max_oxygen_atoms != -1: if struct.get_num_atoms('O') > max_oxygen_atoms: - return (True, f"Exceeded maximumOxygenAtoms: {struct.get_num_atoms('O')} > {max_oxygen_atoms}") + return f"Exceeded maximumOxygenAtoms: {struct.get_num_atoms('O')} > {max_oxygen_atoms}" max_nitrogen_atoms = species_constraints.get('maximumNitrogenAtoms', -1) if max_nitrogen_atoms != -1: if struct.get_num_atoms('N') > max_nitrogen_atoms: - return (True, f"Exceeded maximumNitrogenAtoms: {struct.get_num_atoms('N')} > {max_nitrogen_atoms}") + return f"Exceeded maximumNitrogenAtoms: {struct.get_num_atoms('N')} > {max_nitrogen_atoms}" max_silicon_atoms = species_constraints.get('maximumSiliconAtoms', -1) if max_silicon_atoms != -1: if struct.get_num_atoms('Si') > max_silicon_atoms: - return (True, f"Exceeded maximumSiliconAtoms: {struct.get_num_atoms('Si')} > {max_silicon_atoms}") + return f"Exceeded maximumSiliconAtoms: {struct.get_num_atoms('Si')} > {max_silicon_atoms}" max_sulfur_atoms = species_constraints.get('maximumSulfurAtoms', -1) if max_sulfur_atoms != -1: if struct.get_num_atoms('S') > max_sulfur_atoms: - return (True, f"Exceeded maximumSulfurAtoms: {struct.get_num_atoms('S')} > {max_sulfur_atoms}") + return f"Exceeded maximumSulfurAtoms: {struct.get_num_atoms('S')} > {max_sulfur_atoms}" max_heavy_atoms = species_constraints.get('maximumHeavyAtoms', -1) if max_heavy_atoms != -1: heavy_atoms = struct.get_num_atoms() - struct.get_num_atoms('H') if heavy_atoms > max_heavy_atoms: - return (True, f"Exceeded maximumHeavyAtoms: {heavy_atoms} > {max_heavy_atoms}") + return f"Exceeded maximumHeavyAtoms: {heavy_atoms} > {max_heavy_atoms}" max_surface_sites = species_constraints.get('maximumSurfaceSites', -1) if max_surface_sites != -1: if struct.get_num_atoms('X') > max_surface_sites: - return (True, f"Exceeded maximumSurfaceSites: {struct.get_num_atoms('X')} > {max_surface_sites}") + return f"Exceeded maximumSurfaceSites: {struct.get_num_atoms('X')} > {max_surface_sites}" max_surface_bond_order = species_constraints.get('maximumSurfaceBondOrder', -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, f"Exceeded maximumSurfaceBondOrder at site: {site.get_total_bond_order()} > {max_surface_bond_order}") + return f"Exceeded maximumSurfaceBondOrder at site: {site.get_total_bond_order()} > {max_surface_bond_order}" max_radicals = species_constraints.get('maximumRadicalElectrons', -1) if max_radicals != -1: if struct.get_radical_count() > max_radicals: - return (True, f"Exceeded maximumRadicalElectrons: {struct.get_radical_count()} > {max_radicals}") + return f"Exceeded maximumRadicalElectrons: {struct.get_radical_count()} > {max_radicals}" max_carbenes = species_constraints.get('maximumSingletCarbenes', 1) - if max_radicals != -1: + if max_carbenes != -1: if struct.get_singlet_carbene_count() > max_carbenes: - return (True, f"Exceeded maximumSingletCarbenes: {struct.get_singlet_carbene_count()} > {max_carbenes}") + return f"Exceeded maximumSingletCarbenes: {struct.get_singlet_carbene_count()} > {max_carbenes}" max_carbene_radicals = species_constraints.get('maximumCarbeneRadicals', 0) if max_carbene_radicals != -1: if struct.get_singlet_carbene_count() > 0 and struct.get_radical_count() > max_carbene_radicals: - return (True, f"Exceeded maximumCarbeneRadicals: {struct.get_radical_count()} > {max_carbene_radicals}") + return f"Exceeded maximumCarbeneRadicals: {struct.get_radical_count()} > {max_carbene_radicals}" - return (False, None) + return False diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 9f47187d037..09553473994 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -1654,8 +1654,8 @@ def _generate_product_structures(self, reactant_structures, maps, forward, relab for struct in product_structures: if self.is_molecule_forbidden(struct): raise ForbiddenStructureException() - failed, reason = fails_species_constraints(struct) - if failed: + if fails_species_constraints(struct): + reason = fails_species_constraints(struct) raise ForbiddenStructureException( "Species constraints forbids product species {0}. Please " "reformulate constraints, or explicitly " diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index b3630846bac..f673912f7fd 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -689,11 +689,11 @@ def initialize(self, **kwargs): "Input species {0} is globally forbidden. You may explicitly " "allow it by adding 'input species' to the `generatedSpeciesConstraints` `allowed` list.".format(spec.label) ) - failed, reason = fails_species_constraints(spec) - if failed: + if fails_species_constraints(spec): if "allowed" in self.species_constraints and "input species" in self.species_constraints["allowed"]: self.species_constraints["explicitlyAllowedMolecules"].append(spec.molecule[0]) else: + reason = fails_species_constraints(spec) raise ForbiddenStructureException( "Species constraints forbids input species {0}. Please " "reformulate constraints, remove the species, or explicitly " diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index feaf92b5eb0..607d4eecae7 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -1673,11 +1673,11 @@ def add_seed_mechanism_to_core(self, seed_mechanism, react=False, requires_rms=F "found in a seed mechanism or reaction " "library.".format(spec.label, seed_mechanism.label) ) - failed, reason = fails_species_constraints(spec) - if failed: + if fails_species_constraints(spec): if "allowed" in rmg.species_constraints and "seed mechanisms" in rmg.species_constraints["allowed"]: rmg.species_constraints["explicitlyAllowedMolecules"].extend(spec.molecule) else: + reason = fails_species_constraints(spec) raise ForbiddenStructureException( "Species constraints forbids species {0} from seed mechanism {1}." " Please reformulate constraints, remove the species, or" @@ -1801,11 +1801,11 @@ def add_reaction_library_to_edge(self, reaction_library, requires_rms=False): "inert unless found in a seed mechanism or reaction " "library.".format(spec.label, reaction_library.label) ) - failed, reason = fails_species_constraints(spec) - if failed: + if fails_species_constraints(spec): if "allowed" in rmg.species_constraints and "reaction libraries" in rmg.species_constraints["allowed"]: rmg.species_constraints["explicitlyAllowedMolecules"].extend(spec.molecule) else: + reason = fails_species_constraints(spec) raise ForbiddenStructureException( "Species constraints forbids species {0} from reaction library " "{1}. Please reformulate constraints, remove the species, or " diff --git a/test/rmgpy/constraintsTest.py b/test/rmgpy/constraintsTest.py index 9db23f33a0c..83097d26d60 100644 --- a/test/rmgpy/constraintsTest.py +++ b/test/rmgpy/constraintsTest.py @@ -83,8 +83,7 @@ def test_constraints_not_loaded(self, mock_logging): rmgpy.rmg.input.rmg = None mol = Molecule(smiles="C") - failed, _ = fails_species_constraints(mol) - assert not failed + assert not fails_species_constraints(mol) mock_logging.debug.assert_called_with("Species constraints could not be found.") @@ -97,80 +96,67 @@ def test_species_input(self): """ spc = Species().from_smiles("C") - failed, _ = fails_species_constraints(spc) - assert not failed + assert not fails_species_constraints(spc) def test_explicitly_allowed_molecules(self): """ Test that we can explicitly allow molecules in species constraints. """ mol = Molecule(smiles="CCCC") - failed, _ = fails_species_constraints(mol) - assert failed + assert fails_species_constraints(mol) self.rmg.species_constraints["explicitlyAllowedMolecules"] = [Molecule(smiles="CCCC")] - failed, _ = fails_species_constraints(mol) - assert not failed + assert not fails_species_constraints(mol) def test_carbon_constraint(self): """ Test that we can constrain the max number of carbon atoms. """ mol1 = Molecule(smiles="CC") - failed, _ = fails_species_constraints(mol1) - assert not failed + assert not fails_species_constraints(mol1) mol2 = Molecule(smiles="CCC") - failed, _ = fails_species_constraints(mol2) - assert failed + assert fails_species_constraints(mol2) def test_oxygen_constraint(self): """ Test that we can constrain the max number of oxygen atoms. """ mol1 = Molecule(smiles="C=O") - failed, _ = fails_species_constraints(mol1) - assert not failed + assert not fails_species_constraints(mol1) mol2 = Molecule(smiles="OC=O") - failed, _ = fails_species_constraints(mol2) - assert failed + assert fails_species_constraints(mol2) def test_nitrogen_constraint(self): """ Test that we can constrain the max number of nitrogen atoms. """ mol1 = Molecule(smiles="CN") - failed, _ = fails_species_constraints(mol1) - assert not failed + assert not fails_species_constraints(mol1) mol2 = Molecule(smiles="NCN") - failed, _ = fails_species_constraints(mol2) - assert failed + assert fails_species_constraints(mol2) def test_silicon_constraint(self): """ Test that we can constrain the max number of silicon atoms. """ mol1 = Molecule(smiles="[SiH4]") - failed, _ = fails_species_constraints(mol1) - assert not failed + assert not fails_species_constraints(mol1) mol2 = Molecule(smiles="[SiH3][SiH3]") - failed, _ = fails_species_constraints(mol2) - assert failed + assert fails_species_constraints(mol2) def test_sulfur_constraint(self): """ Test that we can constrain the max number of sulfur atoms. """ mol1 = Molecule(smiles="CS") - failed, _ = fails_species_constraints(mol1) - assert not failed + assert not fails_species_constraints(mol1) mol2 = Molecule(smiles="SCS") - failed, _ = fails_species_constraints(mol2) - assert failed + assert fails_species_constraints(mol2) def test_surface_site_constraint(self): """ @@ -222,17 +208,13 @@ def test_surface_site_constraint(self): self.rmg.species_constraints["maximumCarbonAtoms"] = 3 self.rmg.species_constraints["maximumHeavyAtoms"] = 6 - failed, _ = fails_species_constraints(mol_1site) - assert not failed + assert not fails_species_constraints(mol_1site) - failed, _ = fails_species_constraints(mol_2site) - assert not failed + assert not fails_species_constraints(mol_2site) - failed, _ = fails_species_constraints(mol_3site_vdW) - assert failed + assert fails_species_constraints(mol_3site_vdW) - failed, _ = fails_species_constraints(mol_3site) - assert failed + assert fails_species_constraints(mol_3site) self.rmg.species_constraints["maximumCarbonAtoms"] = max_carbon self.rmg.species_constraints["maximumHeavyAtoms"] = max_heavy_atoms @@ -247,32 +229,27 @@ def test_surface_bond_order_constraint(self): 2 X u0 p0 c0 {1,Q} """ ) - failed, _ = fails_species_constraints(mol_1site) - assert failed + assert fails_species_constraints(mol_1site) def test_heavy_constraint(self): """ Test that we can constrain the max number of heavy atoms. """ mol1 = Molecule(smiles="CCO") - failed, _ = fails_species_constraints(mol1) - assert not failed + assert not fails_species_constraints(mol1) mol2 = Molecule(smiles="CCN=O") - failed, _ = fails_species_constraints(mol2) - assert failed + assert fails_species_constraints(mol2) def test_radical_constraint(self): """ Test that we can constrain the max number of radical electrons. """ mol1 = Molecule(smiles="[CH2][CH2]") - failed, _ = fails_species_constraints(mol1) - assert not failed + assert not fails_species_constraints(mol1) mol2 = Molecule(smiles="[CH2][CH][CH2]") - failed, _ = fails_species_constraints(mol2) - assert failed + assert fails_species_constraints(mol2) def test_carbene_constraint(self): """ @@ -285,8 +262,7 @@ def test_carbene_constraint(self): 3 H u0 p0 c0 {1,S} """ ) - failed, _ = fails_species_constraints(mol1) - assert not failed + assert not fails_species_constraints(mol1) mol2 = Molecule().from_adjacency_list( """ @@ -296,8 +272,7 @@ def test_carbene_constraint(self): 4 H u0 p0 c0 {3,S} """ ) - failed, _ = fails_species_constraints(mol2) - assert failed + assert fails_species_constraints(mol2) def test_carbene_radical_constraint(self): """ @@ -310,8 +285,7 @@ def test_carbene_radical_constraint(self): 3 H u0 p0 c0 {1,S} """ ) - failed, _ = fails_species_constraints(mol1) - assert not failed + assert not fails_species_constraints(mol1) mol2 = Molecule().from_adjacency_list( """ @@ -322,5 +296,4 @@ def test_carbene_radical_constraint(self): 5 H u0 p0 c0 {3,S} """ ) - failed, _ = fails_species_constraints(mol2) - assert failed + assert fails_species_constraints(mol2)