Atom Mapping Fixup (#672)

This PR solves a problem with the atom mapping algorithm in ARC, in
which multiple scissions are required on the same species. A new method
is introduced and tested. In addition, the make_bond_change function can
now add a double bond between atoms that one of them does not have
radicals. Additional fix ups are introduced.

arc/mapping/driver.py

@@ -243,10 +243,11 @@ def map_rxn(rxn: 'ARCReaction',
     #step 3:
     label_species_atoms(reactants)
     label_species_atoms(products)
-    r_cuts, p_cuts = cut_species_for_mapping(reactants, products,loc_r,loc_p)
+    r_cuts = cut_species_for_mapping(reactants, loc_r)
+    p_cuts = cut_species_for_mapping(products, loc_p)
 
     try:
-        make_bond_changes(rxn,r_cuts,r_label_dict)
+        make_bond_changes(rxn, r_cuts, r_label_dict)
     except (ValueError, IndexError, ActionError):
         return None
 

arc/mapping/driver_test.py

@@ -21,7 +21,7 @@
 from itertools import permutations
 
 
-class TestMapping(unittest.TestCase):
+class TestMappingDriver(unittest.TestCase):
     """
     Contains unit tests for the mapping module.
     """

arc/mapping/engine.py

@@ -1055,13 +1055,12 @@ def prepare_reactants_and_products_for_scissors(rxn: 'ARCReaction',
                 else:
                     loc_r[location] += 1
                     reactants[location] = ARCSpecies(label="".join(sorted(
-                        [key_by_val(r_label_dict, r_label_dict[broken_bond[1]]),
-                         key_by_val(p_label_dict, p_label_dict[broken_bond[3]])])),
-                    mol = reactant.mol.copy(deep=True),
-                    xyz = reactant.get_xyz(),
-                    bdes = [(r_label_dict[broken_bond[1]] + 1 - index,
-                             r_label_dict[broken_bond[3]] + 1 - index)])
-
+                                                        [key_by_val(r_label_dict, r_label_dict[broken_bond[1]]),
+                                                        key_by_val(p_label_dict, p_label_dict[broken_bond[3]])])),
+                                                     mol = reactant.mol.copy(deep=True),
+                                                     bdes = [(r_label_dict[broken_bond[1]] + 1 - index,
+                                                              r_label_dict[broken_bond[3]] + 1 - index)])
+                    reactants[location].final_xyz = reactant.get_xyz()
                     for mol1, mol2 in zip(reactants[location].mol.atoms, reactant.mol.atoms):
                         mol1.label = mol2.label
                     break
@@ -1103,32 +1102,42 @@ def prepare_reactants_and_products_for_scissors(rxn: 'ARCReaction',
 def make_bond_changes(rxn: 'ARCReaction',
                       r_cuts: List[ARCSpecies],
                       r_label_dict: Dict):
-
     """
     Makes the bond change before matching the reactants and products
 
     Ags:
         rxn: An ARCReaction object
         r_cuts: the cut products
         r_label_dict: the dictionary object the find the relevant location.
-
-    Returns:
-        None if there are no bond changes, r_cuts with a muted ARCSpecies in the relevant location
     """
+
     for action in rxn.family.forward_recipe.actions:
         if action[0].lower() == "CHANGE_BOND".lower():
             indicies = r_label_dict[action[1]],r_label_dict[action[3]]
             for r_cut in r_cuts:
                 if indicies[0] in [int(atom.label) for atom in r_cut.mol.atoms] and indicies[1] in [int(atom.label) for atom in r_cut.mol.atoms]:
-                    atom1, atom2 = 0,0
+                    atom1, atom2 = 0, 0
                     for atom in r_cut.mol.atoms:
                         if int(atom.label) == indicies[0]:
                             atom1 = atom
                         elif int(atom.label) == indicies[1]:
                             atom2 = atom
-                    atom1.decrement_radical()
-                    atom2.decrement_radical()
-                    r_cut.mol.get_bond(atom1,atom2).order += action[2]
+                    if atom1.radical_electrons == 0:
+                        atom1.lone_pairs -= 1
+                        atom2.lone_pairs += 1
+                        atom1.charge += 1
+                        atom2.charge -= 1
+                        atom2.radical_electrons -= 2
+                    elif atom2.radical_electrons == 0:
+                        atom2.lone_pairs -= 1
+                        atom1.lone_pairs += 1
+                        atom2.charge += 1
+                        atom1.charge -= 1
+                        atom1.radical_electrons -= 2
+                    else:    
+                        atom1.decrement_radical()
+                        atom2.decrement_radical()
+                    r_cut.mol.get_bond(atom1, atom2).order += action[2]
                     r_cut.mol.update()
 
 
@@ -1154,83 +1163,66 @@ def assign_labels_to_products(rxn: 'ARCReaction',
             atom_index+=1
 
 
-def cut_species_for_mapping(reactants: List[ARCSpecies],
-                            products: List[ARCSpecies],
-                            loc_r: List[int],
-                            loc_p: List[int],
-                            ) -> Optional[Tuple[List[ARCSpecies], List[ARCSpecies]]]:
+def multiple_cut_on_species(spc, bdes):
     """
-    A function for scissoring the reactants and products, as a preparation for atom mapping.
-
+    A function that recursively calles an is called by cut_species_for_mapping.
+    Is used to cut a species and reassign the other BDE's to the other cuts, and recalling cut_species_for_mapping.
     Args:
-        reactants: A list of the ARCSpecies for scission
-        products: A list of the ARCSpecies for scission
-        loc_r: A list of the location and number of cuts that is required.
-        loc_p: A list of the location and number of cuts that is required.
+        spc (ARCSpecies): a species with more then one BDE's (marks for scission).
+        bdes (list(tuple(int))): the required BDE's.
+    Returns:
+        list(ARCSpecies): a list of the cut products.
+    """
+    bdes = spc.bdes
+    spc.bdes = [bdes[0]]
+    bdes = bdes[1:]
+    try:
+        cuts = spc.scissors()
+    except SpeciesError as e:
+        return None
+    for species in cuts:
+        species.final_xyz = species.get_xyz(generate=False)
+        indinces = [int(atom.label) for atom in species.mol.copy(deep=True).atoms]
+        new_bdes = list()
+        for bde in bdes:
+            if bde[0]-1 in indinces and bde[1]-1 in indinces:
+                new_bdes.append((indinces.index(bde[0]-1) + 1, indinces.index(bde[1]-1) + 1))
+        species.bdes = new_bdes
+    return cut_species_for_mapping(cuts, [len(species.bdes or list()) for species in cuts])
+
 
+def cut_species_for_mapping(species, locs):
+    """
+    A function for performing the necessary scission of species for mapping purposes. Can perform appropriate scission of multiple bonds at once.
+    Args:
+        species (list(ARCSpecies)): the species (reactants or products), marked for scission.
+        locs (list(int)): the number of cuts that is required for each species.
     Returns:
-        A list of scissored reactants and products.
+        list(ARCSpecies): a list of the cut products.
     """
-    r_cuts, p_cuts=list(), list()
-    for index, reactant in zip(loc_r, reactants):
-        if index==1:
+    cuts = list()
+    for spc, loc in zip(species, locs):
+        spc.final_xyz = spc.get_xyz()
+        if spc.mol.copy(deep=True).smiles == "[H][H]" and loc != 0: # scissors should return one species
+            labels = [atom.label for atom in spc.mol.copy(deep=True).atoms]
             try:
-                reactant.final_xyz = reactant.get_xyz()
-                cuts = reactant.scissors()
-                r_cuts += cuts
+                H1 = spc.scissors()[0]
             except SpeciesError:
                 return None
-        elif index>1:
-            bdes = reactant.bdes
-            new_r = ARCSpecies(label="scissors", mol=reactant.mol.copy(deep=True))
-            for bde in bdes:
-                new_r.bdes = [bde]
-                new_r.final_xyz = new_r.get_xyz()
-                try:
-                    cuts=new_r.scissors()
-                except SpeciesError:
-                    return None
-                if len(cuts) == 1:
-                    new_r = cuts[0]
-                else:
-                    new_r, second = find_main_cut_product(cuts, reactant, bde)
-                    r_cuts += [second]
-            r_cuts += [new_r]
-        else:
-            r_cuts.append(reactant)
-
-    for index, product in zip(loc_p, products):
-        if index==1:
+            H2 = H1.copy()
+            H2.mol.atoms[0].label = labels[0] if H1.mol.atoms[0].label != labels[0] else labels[1]
+            cuts += [H1, H2]
+        if loc == 0:
+            cuts += [spc]
+        elif loc == 1:
             try:
-                product.final_xyz = product.get_xyz()
-                cuts = product.scissors()
-                if len(cuts) == 1: #only H2 and cyclic species for now, todo: modify to include cyclic species.
-                    cuts.append(ARCSpecies(label= cuts[0].label, mol=cuts[0].mol.copy(deep=True)))
-                    labels = [atom.label for atom in product.mol.atoms]
-                    cuts[-1].mol.atoms[0].label = labels[1] if cuts[0].mol.atoms[0].label == labels[0] else labels[0]
-                p_cuts += cuts
+                cuts += spc.scissors()
             except SpeciesError:
                 return None
-        elif index > 1:
-            bdes = product.bdes
-            new_p = ARCSpecies(label="scissors", mol=product.mol.copy(deep=True))
-            for bde in bdes:
-                new_p.bdes = [bde]
-                new_p.final_xyz = new_p.get_xyz()
-                try:
-                    cuts = new_p.scissors()
-                except SpeciesError:
-                    return None
-                if len(cuts) == 1:
-                    new_p = cuts[0]
-                else:
-                    new_p, second = find_main_cut_product(cuts, product, bde)
-                    p_cuts += [second]
-            p_cuts += [new_p]
         else:
-            p_cuts.append(product)
-
-    return r_cuts, p_cuts
+            bdes = spc.bdes
+            cuts += multiple_cut_on_species(spc, bdes)
+    return cuts
 
 
 def find_main_cut_product(cuts: List["ARCSpecies"],