Success in YOBCUO

Author: choglass

cell2mol/cell_reconstruction.py

@@ -2,8 +2,11 @@
 import itertools
 from cell2mol.cell_operations import translate
 from cell2mol.other import additem, absolute_value
-from cell2mol.connectivity import compare_species, count_species, split_species
+from cell2mol.connectivity import compare_species, count_species, split_species, arrange_data_for_reorder
+from cell2mol.hungarian import reorder
 from cell2mol.elementdata import ElementData
+from cell2mol.read_write import writexyz
+
 elemdatabase = ElementData()
 
 #######################################################
@@ -208,6 +211,9 @@ def fragments_reconstruct(moleclist: list, fraglist: list, Hlist: list, refmolec
         moleclist.extend(finalmols)
         print(f"{finalmols=}")
         print(f"{remfrag=}")
+        for i, g in enumerate(finalmols):
+            if debug == 1: writexyz("/Users/ycho/cell2mol/cell2mol/test/YOBCUO/", f"reorder_molec_{i}.xyz", g.labels, g.coord)
+
         if len(remfrag) > 0: Warning = True;  print("FRAG_RECONSTRUCT. Remaining after Hydrogen reconstruction",remfrag)
         else:                Warning = False; print("FRAG_RECONSTRUCT. No remaining Molecules after Hydrogen reconstruction")
     elif len(remfrag) > 0 and len(Hlist) == 0:
@@ -221,6 +227,7 @@ def fragments_reconstruct(moleclist: list, fraglist: list, Hlist: list, refmolec
     # The former were identified in the cell.get_moleclist() function, and have cell as parent. 
     # The latter have been constructed by merging fragments, and do not have cell as parent, but have the cell_indices stored in mol.cell_indices
     # Here we homogenize the situation by adding the cell_indices variable to all molecules
+    # TODO : why cell_indices is needed?
     for mol in moleclist:
         if not hasattr(mol,"cell_indices"): 
             if mol.check_parent("cell"):
@@ -466,6 +473,8 @@ def sequential(fragmentlist: list, refmoleclist: list, cellvec: list, factor: fl
 
 #######################################################
 def combine(tobemerged: list, references: list, cellvec: list, threshold_tmat: float, cov_factor: float, metal_factor: float, debug: int=0):
+    from cell2mol.classes import molecule
+
     goodlist = []   ## List of molecules coming from the two fragments received
     avglist = []    ## List of bigger fragments coming from the two fragments received
     badlist = []    ## List of fragments as they entered the function
@@ -488,12 +497,33 @@ def combine(tobemerged: list, references: list, cellvec: list, threshold_tmat: f
         found = False 
         for ref in references:
             if not found: 
-                issame = compare_species(newmolec, ref)
-                if issame:    ## Then is a molecule that appears in the reference list 
-                    found = True 
-                    newmolec.subtype = ref.subtype
-                    goodlist.append(newmolec)
-                    if debug >= 1: print(f"COMBINE: Fragment {newmolec.formula} added to goodlist with {newmolec.cell_indices=}")
+                if (newmolec.natoms == ref.natoms) and (newmolec.eleccount == ref.eleccount) and (newmolec.formula == ref.formula):
+                    dummy1, dummy2, map12 = reorder(ref.labels, newmolec.labels, ref.coord, newmolec.coord)
+                    
+                    reordered_labels = [newmolec.labels[i] for i in map12]
+                    reordered_coord  = [newmolec.coord[i] for i in map12]
+                    reordered_radii  = [newmolec.radii[i] for i in map12]
+                    reordered_frac_cood = [newmolec.frac_coord[i] for i in map12]
+                    reordered_cell_indices = [newmolec.cell_indices[i] for i in map12]
+
+                    reordered_newmolec = molecule(reordered_labels, reordered_coord, reordered_radii)
+                    reordered_newmolec.cell_indices = reordered_cell_indices
+                    reordered_newmolec.set_fractional_coord(reordered_frac_cood)
+                    reordered_newmolec.set_atoms(create_adjacencies=True, debug=2)
+                    if reordered_newmolec.iscomplex: 
+                        reordered_newmolec.split_complex()
+                        reordered_newmolec.get_hapticity(debug=debug)
+                        for lig in reordered_newmolec.ligands:
+                            lig.get_denticity(debug=debug)
+
+                    print(f"{reordered_newmolec=}")
+
+                    issame = compare_species(reordered_newmolec, ref, debug=2)
+                    if issame:    ## Then is a molecule that appears in the reference list 
+                        found = True 
+                        reordered_newmolec.subtype = ref.subtype
+                        goodlist.append(reordered_newmolec)
+                        if debug >= 1: print(f"COMBINE: Fragment {reordered_newmolec.formula} added to goodlist with {reordered_newmolec.cell_indices=}")
         if not found:        ## Then it is a fragment. A bigger one, but still a fragment
             newmolec.subtype = "Rec. Fragment"
             avglist.append(newmolec)

cell2mol/charge_assignment.py

@@ -864,40 +864,6 @@ def prepare_unresolved(unique_indices: list, unique_species: list, distributions
 
     return list_molecules, list_indices, list_options
 
-#######################################################
-def arrange_data_for_reorder(reference: object, target: object, debug: int=0):
-    # To do the reorder, we create new tags that include as much information as possible.
-    # Ideally, we aim to include the label + the connectivity + the metal connectivity
-    t_totconnec = 0
-    t_totmconnec = 0
-    for a in target.atoms:
-        t_totconnec  += a.connec
-        t_totmconnec += a.mconnec
-    r_totconnec = 0
-    r_totmconnec = 0
-    for a in reference.atoms:
-        r_totconnec  += a.connec
-        r_totmconnec += a.mconnec
-    if t_totconnec == r_totconnec:   useconec = True
-    else:                            useconec = False
-    if t_totmconnec == r_totmconnec: usemconec = True
-    else:                            usemconec = False
-    # For target
-    target_data = []
-    for a in target.atoms:
-        data = a.label
-        if useconec:  data += str(a.connec)
-        if usemconec: data += str(a.mconnec)
-        target_data.append(data)
-    # For reference
-    ref_data = []
-    for a in reference.atoms:
-        data = a.label
-        if useconec:  data += str(a.connec)
-        if usemconec: data += str(a.mconnec)
-        ref_data.append(data)
-    return ref_data, target_data
-
 #######################################################    
 def set_charges_create_bonds (specie, unique_indices, unique_species, final_charge_distribution, debug):
 

cell2mol/classes.py

@@ -1356,7 +1356,7 @@ def reconstruct(self, cov_factor: float=None, metal_factor: float=None, debug: i
                 newmolec = molecule(mol.labels, mol.coord)
                 newmolec.origin = "cell.reconstruct"
                 newmolec.set_atoms(create_adjacencies=True, debug=debug)
-                newmolec.add_parent(self,mol.cell_indices) 
+                newmolec.add_parent(self, mol.cell_indices) 
                 newmolec.set_fractional_coord(mol.frac_coord)
                 if newmolec.iscomplex: newmolec.split_complex()
                 self.moleclist.append(newmolec)         

cell2mol/connectivity.py

@@ -6,6 +6,7 @@
 from typing import Tuple
 from cell2mol.other import inv
 from cell2mol.elementdata import ElementData
+from cell2mol.read_write import writexyz
 elemdatabase = ElementData()
 
 #######################################################
@@ -374,6 +375,12 @@ def compare_species(mol1, mol2, check_coordinates: bool=False, debug: int=0):
         if debug == 1: print(mol2.formula)
         if debug == 2: print(mol1)
         if debug == 2: print(mol2)
+        if debug == 2: print(mol1.labels)
+        if debug == 2: print(mol2.labels)
+
+    if debug == 2: writexyz("/Users/ycho/cell2mol/cell2mol/test/YOBCUO/", "reorder_molec.xyz", mol1.labels, mol1.coord)
+    if debug == 2: writexyz("/Users/ycho/cell2mol/cell2mol/test/YOBCUO/", "ref.xyz", mol2.labels, mol2.coord)
+
     # a pair of species is compared on the basis of:
     # 1) the total number of atoms
     if (mol1.natoms != mol2.natoms): 
@@ -396,10 +403,16 @@ def compare_species(mol1, mol2, check_coordinates: bool=False, debug: int=0):
     # 4) the number of adjacencies between each pair of element types
     if not hasattr(mol1,"adj_types"):     mol1.set_adj_types()
     if not hasattr(mol2,"adj_types"):     mol2.set_adj_types()
+    if debug == 2: print(f"{mol1.adj_types=}")
+    if debug == 2: print(f"{mol2.adj_types=}")
+    if debug == 2: np.save("/Users/ycho/cell2mol/cell2mol/test/YOBCUO/adj_types1.npy", mol1.adj_types)
+    if debug == 2: np.save("/Users/ycho/cell2mol/cell2mol/test/YOBCUO/adj_types2.npy", mol2.adj_types)
     for kdx, elem in enumerate(mol1.adj_types):
         for ldx, elem2 in enumerate(elem):
             if elem2 != mol2.adj_types[kdx, ldx]: 
                 if debug > 0: print(f"COMPARE_SPECIES. FALSE, different adjacency count")
+                if debug > 0: print(f"{kdx=} {ldx=} {elem=} {elem2=} {mol2.adj_types[kdx]=} {mol2.adj_types[kdx, ldx]=}")
+                if debug > 0: print(f"{mol1.labels[kdx]=} {mol1.labels[ldx]=} {mol2.labels[kdx]=} {mol2.labels[ldx]=}")
                 return False
 
     if check_coordinates:
@@ -411,4 +424,35 @@ def compare_species(mol1, mol2, check_coordinates: bool=False, debug: int=0):
     return True
 
 #################################
-
+def arrange_data_for_reorder(reference: object, target: object, debug: int=0):
+    # To do the reorder, we create new tags that include as much information as possible.
+    # Ideally, we aim to include the label + the connectivity + the metal connectivity
+    t_totconnec = 0
+    t_totmconnec = 0
+    for a in target.atoms:
+        t_totconnec  += a.connec
+        t_totmconnec += a.mconnec
+    r_totconnec = 0
+    r_totmconnec = 0
+    for a in reference.atoms:
+        r_totconnec  += a.connec
+        r_totmconnec += a.mconnec
+    if t_totconnec == r_totconnec:   useconec = True
+    else:                            useconec = False
+    if t_totmconnec == r_totmconnec: usemconec = True
+    else:                            usemconec = False
+    # For target
+    target_data = []
+    for a in target.atoms:
+        data = a.label
+        if useconec:  data += str(a.connec)
+        if usemconec: data += str(a.mconnec)
+        target_data.append(data)
+    # For reference
+    ref_data = []
+    for a in reference.atoms:
+        data = a.label
+        if useconec:  data += str(a.connec)
+        if usemconec: data += str(a.mconnec)
+        ref_data.append(data)
+    return ref_data, target_data