Pdep spin conservation

rmgpy/data/kinetics/database.py

@@ -436,6 +436,7 @@ def generate_reactions(self, reactants, products=None, only_families=None, reson
         reaction_list = []
         if only_families is None:
             reaction_list.extend(self.generate_reactions_from_libraries(reactants, products))
+
         reaction_list.extend(self.generate_reactions_from_families(reactants, products,
                                                                    only_families=only_families, resonance=resonance))
         return reaction_list

rmgpy/data/kinetics/family.py

@@ -213,6 +213,33 @@ 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_spin_string =  "+".join(sorted([str(reactant.multiplicity) for reactant in self.reactants]))
+            reactant_combined_spin = allowed_spin[reactant_spin_string]
+        else:
+            return None
+        if len(self.products) == 1:
+            product_combined_spin = {self.products[0].multiplicity}
+        elif len(self.products) == 2:
+            product_spin_string = "+".join(sorted([str(product.multiplicity) for product in self.products]))
+            product_combined_spin = allowed_spin[product_spin_string]
+        else:
+            return None
+        return reactant_combined_spin, product_combined_spin
 
     def apply_solvent_correction(self, solvent):
         """
@@ -1653,6 +1680,7 @@ def _generate_product_structures(self, reactant_structures, maps, forward, relab
         # Apply the generated species constraints (if given)
         for struct in product_structures:
             if self.is_molecule_forbidden(struct):
+                # logging.info(f"{str(struct)} is forbidden!")
                 raise ForbiddenStructureException()
             if fails_species_constraints(struct):
                 raise ForbiddenStructureException()
@@ -1678,7 +1706,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`
@@ -1713,7 +1741,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.info("Did not create the following reaction, which violates conservation of spin...")
+                logging.info(str(reaction))
+                return None
         return reaction
 
     def _match_reactant_to_template(self, reactant, template_reactant):
@@ -1963,7 +1995,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
@@ -1982,6 +2014,7 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson
             reactant_num = self.product_num
 
         if self.auto_generated and reactant_num != len(reactants):
+            # logging.info("self.auto_generated and reactant_num != len(reactants)")
             return []
 
         if len(reactants) > len(template.reactants):
@@ -2019,7 +2052,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
@@ -2062,7 +2097,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)
 
@@ -2086,8 +2123,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)
 
@@ -2140,7 +2178,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:
@@ -2205,7 +2245,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)
 
@@ -4897,3 +4939,19 @@ def get_site_solute_data(rxn):
         return site_data
     else:
         return None
+
+allowed_spin_violation_families =['1,2-Birad_to_alkene','1,4_Cyclic_birad_scission','1,4_Linear_birad_scission']
+allowed_spin = {
+    "1+1": set([1]),
+    "1+2": set([2]),
+    "1+3": set([3]),
+    "1+4": set([4]),
+    "1+5": set([5]),
+    "2+2": set([1,3]),
+    "2+3": set([2,4]),
+    "2+4": set([3,5]),
+    "2+5": set([4,6]),
+    "3+3": set([1,3,5]),
+    "3+4": set([2,4,6]),
+    "3+5": set([3,5,7]),
+}

rmgpy/molecule/molecule.py

@@ -59,6 +59,7 @@
 from rmgpy.molecule.kekulize import kekulize
 from rmgpy.molecule.pathfinder import find_shortest_path
 from rmgpy.molecule.fragment import CuttingLabel
+from rmgpy.molecule.util import generate_closed_shell_singlet, generate_singlet_diradicals
 
 ################################################################################
 
@@ -3005,6 +3006,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