Pdep spin conservation

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):
@@ -1807,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 
@@ -1863,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
@@ -1906,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)
 
@@ -1930,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)
 
@@ -1984,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:
@@ -2049,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)
 
@@ -4495,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']

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

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

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