resolved merge conflict

Author: donerancl

rmgpy/data/kinetics/family.py

@@ -213,6 +213,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
 
     def apply_solvent_correction(self, solvent):
         """
@@ -1678,7 +1703,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 +1738,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):

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