From 49e67839cb217c00e3c44dae16238cab8e1fa694 Mon Sep 17 00:00:00 2001 From: Anna Doner Date: Mon, 12 Aug 2024 16:59:46 -0400 Subject: [PATCH 1/3] reset to main, then added check for reversibility before treating a reaction as pdep --- rmgpy/rmg/model.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index 1ab45dc67c..994bdf5558 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -849,6 +849,8 @@ def process_new_reactions(self, new_reactions, new_species, pdep_network=None, g pdep = rxn.generate_high_p_limit_kinetics() elif any([any([x.is_subgraph_isomorphic(q) for q in self.unrealgroups]) for y in rxn.reactants + rxn.products for x in y.molecule]): pdep = False + elif not rxn.reversible: + pdep = False # If pressure dependence is on, we only add reactions that are not unimolecular; # unimolecular reactions will be added after processing the associated networks From cdbca51754f643fb8b04c67835985fdead8dd210 Mon Sep 17 00:00:00 2001 From: Anna Doner Date: Fri, 23 Aug 2024 14:09:33 -0400 Subject: [PATCH 2/3] add tools to check if reaction template conserves spin, check spin before creating reaction --- rmgpy/data/kinetics/family.py | 33 +++++++++++++++++++++++++++-- rmgpy/molecule/molecule.py | 15 ++++++++++++++ rmgpy/molecule/util.py | 39 +++++++++++++++++++++++++++++++++++ 3 files changed, 85 insertions(+), 2 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index f8818c78ac..e3a85ed645 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -204,6 +204,31 @@ def copy(self): other.is_forward = self.is_forward return other + + def check_if_spin_allowed(self): + # get the combined spin for reactants and products + reactants_combined_spin, products_combined_spin = self.calculate_combined_spin() + # check if there are any matches for combined spin between reactants and products + if reactants_combined_spin.intersection(products_combined_spin) != set([]): + return True + else: + logging.debug(f"Reactants combined spin is {reactants_combined_spin}, but the products combined spin is {products_combined_spin}") + return False + + def calculate_combined_spin(self): + if len(self.reactants) == 1: + reactant_combined_spin = {self.reactants[0].multiplicity} + elif len(self.reactants) == 2: + reactant_combined_spin = set([self.reactants[0].multiplicity + self.reactants[1].multiplicity -1, np.abs(self.reactants[0].multiplicity-self.reactants[1].multiplicity)+1]) + else: + return None + if len(self.products) == 1: + product_combined_spin = {self.products[0].multiplicity} + elif len(self.products) == 2: + product_combined_spin = set([self.products[0].multiplicity + self.products[1].multiplicity -1, np.abs(self.products[0].multiplicity-self.products[1].multiplicity)+1]) + else: + return None + return reactant_combined_spin, product_combined_spin ################################################################################ @@ -1539,7 +1564,7 @@ def is_molecule_forbidden(self, molecule): return False - def _create_reaction(self, reactants, products, is_forward): + def _create_reaction(self, reactants, products, is_forward, check_spin = True): """ Create and return a new :class:`Reaction` object containing the provided `reactants` and `products` as lists of :class:`Molecule` @@ -1565,7 +1590,11 @@ def _create_reaction(self, reactants, products, is_forward): for key, species_list in zip(['reactants', 'products'], [reaction.reactants, reaction.products]): for species in species_list: reaction.labeled_atoms[key] = dict(reaction.labeled_atoms[key], **species.get_all_labeled_atoms()) - + if check_spin: + if not reaction.check_if_spin_allowed(): + logging.warn("Did not create the following reaction, which violates conservation of spin...") + logging.warn(str(reaction)) + return None return reaction def _match_reactant_to_template(self, reactant, template_reactant): diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 03727be792..5522d1dae2 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -58,6 +58,7 @@ from rmgpy.molecule.graph import Vertex, Edge, Graph, get_vertex_connectivity_value from rmgpy.molecule.kekulize import kekulize from rmgpy.molecule.pathfinder import find_shortest_path +from rmgpy.molecule.util import generate_closed_shell_singlet, generate_singlet_diradicals ################################################################################ @@ -2925,6 +2926,20 @@ def get_desorbed_molecules(self): logging.debug("After removing from surface:\n" + desorbed_molecule.to_adjacency_list()) return desorbed_molecules + def update_to_closed_shell_singlet(self): + for i in range(len(self.atoms)): + self.atoms[i].id = i + assert self.multiplicity == 1 and self.get_radical_count()>0 + radical_center_ids = [x.id for x in self.atoms if x.radical_electrons > 0] + # remove radicals from radical centers (2) + for radical_center_id in radical_center_ids: + self.atoms[radical_center_id].decrement_radical() + # add removed radicals (2) to one of the radical sites as a lone pair (1) + self.atoms[radical_center_ids[0]].increment_lone_pairs() + # pick the best resonance structure + self.generate_resonance_structures() + self.reactive = True + # this variable is used to name atom IDs so that there are as few conflicts by # using the entire space of integer objects diff --git a/rmgpy/molecule/util.py b/rmgpy/molecule/util.py index d1093f4582..1c3da1653b 100644 --- a/rmgpy/molecule/util.py +++ b/rmgpy/molecule/util.py @@ -179,3 +179,42 @@ def swap(to_be_swapped, sample): new_partner = (to_be_swapped - sample).pop() return central, original, new_partner + +def generate_closed_shell_singlet(m: Molecule): + for i in range(len(m.atoms)): + m.atoms[i].id = i + assert m.multiplicity == 1 and m.get_radical_count()>0 + radical_center_ids = [x.id for x in m.atoms if x.radical_electrons > 0] + # remove radicals from radical centers (2) + for radical_center_id in radical_center_ids: + m.atoms[radical_center_id].decrement_radical() + # add removed radicals (2) to one of the radical sites as a lone pair (1) + m.atoms[radical_center_ids[0]].increment_lone_pairs() + # pick the best resonance structure + print([x.smiles for x in m.generate_resonance_structures()]) + return m.generate_resonance_structures()[0] + +def generate_singlet_diradicals(m: Molecule): + for i in range(len(m.atoms)): + m.atoms[i].id = i + + singlet_diradicals = [] + for edge in m.get_all_edges(): + + M = m.copy(deep=True) + if edge.get_order_num() > 1: # find a pi bond + atom1_id = edge.atom1.id + atom2_id = edge.atom2.id + M.atoms[atom1_id].increment_radical() # add a radical to each atom of the pi bond + M.atoms[atom2_id].increment_radical() + M.get_bond(M.atoms[atom1_id], M.atoms[atom2_id]).decrement_order() # remove 1 pi bond + potential_singlet_diradicals = M.generate_resonance_structures() # generate resonance structures + + for potential_singlet_diradical in potential_singlet_diradicals: # find all resonance structures with non-neighboring radical sites + radical_center_ids = sorted([x.id for x in potential_singlet_diradical.atoms if x.radical_electrons==1]) + potential_singlet_diradical_edges = potential_singlet_diradical.get_all_edges() + potential_singlet_diradical_edge_ids = [sorted([x.atom1.id, x.atom2.id]) for x in potential_singlet_diradical_edges] + if radical_center_ids not in potential_singlet_diradical_edge_ids: + if potential_singlet_diradical not in singlet_diradicals: + singlet_diradicals.append(potential_singlet_diradical) + return singlet_diradicals \ No newline at end of file From e956eec62b0ea0782e5a5ed512dc9aa355268d1a Mon Sep 17 00:00:00 2001 From: Anna Doner Date: Fri, 23 Aug 2024 17:18:23 -0400 Subject: [PATCH 3/3] added exception families --- rmgpy/data/kinetics/family.py | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index e3a85ed645..22c10f225a 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -1836,7 +1836,7 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson specified reactants and products within this family. Degenerate reactions are returned as separate reactions. """ - + check_spin = True rxn_list = [] # Wrap each reactant in a list if not already done (this is done to @@ -1892,7 +1892,9 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson pass else: if product_structures is not None: - rxn = self._create_reaction(reactant_structures, product_structures, forward) + if self.label in allowed_spin_violation_families: + check_spin = False + rxn = self._create_reaction(reactant_structures, product_structures, forward, check_spin = check_spin) if rxn: rxn_list.append(rxn) # Bimolecular reactants: A + B --> products @@ -1935,7 +1937,9 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson pass else: if product_structures is not None: - rxn = self._create_reaction(reactant_structures, product_structures, forward) + if self.label in allowed_spin_violation_families: + check_spin = False + rxn = self._create_reaction(reactant_structures, product_structures, forward, check_spin = check_spin) if rxn: rxn_list.append(rxn) @@ -1959,8 +1963,9 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson pass else: if product_structures is not None: - rxn = self._create_reaction(reactant_structures, product_structures, - forward) + if self.label in allowed_spin_violation_families: + check_spin = False + rxn = self._create_reaction(reactant_structures, product_structures, forward, check_spin = check_spin) if rxn: rxn_list.append(rxn) @@ -2013,7 +2018,9 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson pass else: if product_structures is not None: - rxn = self._create_reaction(reactant_structures, product_structures, forward) + if self.label in allowed_spin_violation_families: + check_spin = False + rxn = self._create_reaction(reactant_structures, product_structures, forward, check_spin = check_spin) if rxn: rxn_list.append(rxn) else: @@ -2078,7 +2085,9 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson pass else: if product_structures is not None: - rxn = self._create_reaction(reactant_structures, product_structures, forward) + if self.label in allowed_spin_violation_families: + check_spin = False + rxn = self._create_reaction(reactant_structures, product_structures, forward, check_spin = check_spin) if rxn: rxn_list.append(rxn) @@ -4524,3 +4533,5 @@ def average_kinetics(kinetics_list): Ea=(Ea * 0.001, "kJ/mol"), ) return averaged_kinetics + +allowed_spin_violation_families =['1,2-Birad_to_alkene','1,4_Cyclic_birad_scission','1,4_Linear_birad_scission'] \ No newline at end of file