spin conservation check Nov 22 2024

Author: donerancl

rmgpy/data/kinetics/database.py

@@ -429,6 +429,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

@@ -219,13 +219,15 @@ 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])  
+            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_combined_spin = set([self.products[0].multiplicity + self.products[1].multiplicity -1, np.abs(self.products[0].multiplicity-self.products[1].multiplicity)+1])  
+            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
@@ -1539,6 +1541,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()
@@ -1592,8 +1595,8 @@ def _create_reaction(self, reactants, products, is_forward, check_spin = True):
                 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))
+                logging.info("Did not create the following reaction, which violates conservation of spin...")
+                logging.info(str(reaction))
                 return None
         return reaction
 
@@ -1855,6 +1858,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):
@@ -4534,4 +4538,18 @@ def average_kinetics(kinetics_list):
         )
     return averaged_kinetics
 
-allowed_spin_violation_families =['1,2-Birad_to_alkene','1,4_Cyclic_birad_scission','1,4_Linear_birad_scission']
+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]),
+}