Debugging and get_reference_molecules

cell2mol/c2m_driver.py

@@ -25,6 +25,7 @@
 
     # Filenames for output and cell object
     cell_fname   = os.path.join(current_dir, "Cell_{}.cell".format(name))
+    ref_cell_fname = os.path.join(current_dir, "Ref_Cell_{}.cell".format(name))
     output_fname = os.path.join(current_dir, "cell2mol.out")
 
     ##### Deals with the parsed arguments for verbosity ######
@@ -68,8 +69,9 @@
     # Loads the reference molecules and checks_missing_H
     # TODO : reconstruct the unit cell without using reference molecules
     # TODO : reconstruct the unit cell using (only reconstruction of) reference molecules and Space group
-    newcell.get_reference_molecules(ref_labels, ref_fracs, debug=debug) 
-
+    newcell.get_reference_molecules(ref_labels, ref_fracs, debug=2) 
+    newcell.assess_errors()
+    newcell.save(ref_cell_fname)
     ######################
     ### CALLS CELL2MOL ###
     ######################

cell2mol/classes.py

@@ -484,12 +484,13 @@ def get_connected_atoms(self, debug: int=0):
 
     #######################################################
     def get_denticity(self, debug: int=0):
-        if not hasattr(self,"groups"):      self.split_ligand(debug=debug)
+        #if not hasattr(self,"groups"):      self.split_ligand(debug=2)
         if debug > 0: print(f"LIGAND.Get_denticity: checking connectivity of ligand {self.formula}")
         if debug > 0: print(f"LIGAND.Get_denticity: initial connectivity is {len(self.connected_idx)}")
-
+        print(f"{self.groups=}")
         self.denticity = 0
         for g in self.groups:
+            if debug > 0: print(f"LIGAND.Get_denticity: checking denticity of group {g}")
             self.denticity += g.get_denticity()      ## A check is also performed at the group level
         return self.denticity 
 
@@ -499,7 +500,8 @@ def split_ligand(self, debug: int=0):
         # Identify Connected and Unconnected atoms (to the metal)
         if not hasattr(self,"connected_idx"): self.get_connected_idx()
         conn_idx = self.connected_idx
-
+        self.set_atoms()  
+        # if not hasattr(self,"atoms"):       self.set_atoms() 
         # Split the "ligand to obtain the groups
         conn_labels  = extract_from_list(conn_idx, self.labels, dimension=1)
         conn_coord   = extract_from_list(conn_idx, self.coord, dimension=1)
@@ -510,15 +512,24 @@ def split_ligand(self, debug: int=0):
         for b in blocklist:
             if debug > 0: print(f"PREPARING BLOCK: {b}")
             gr_indices = extract_from_list(b, conn_idx, dimension=1)
-            gr_labels  = extract_from_list(b, self.labels, dimension=1)
-            gr_coord   = extract_from_list(b, self.coord, dimension=1)
-            gr_radii   = extract_from_list(b, self.radii, dimension=1)
-            gr_atoms   = extract_from_list(b, self.atoms, dimension=1)
-            newgroup   = group(gr_labels, gr_coord, parent_indices=gr_indices, radii=gr_radii, parent=self)
+            gr_labels  = extract_from_list(gr_indices, self.labels, dimension=1)
+            gr_coord   = extract_from_list(gr_indices, self.coord, dimension=1)
+            gr_radii   = extract_from_list(gr_indices, self.radii, dimension=1)
+            gr_atoms   = extract_from_list(gr_indices, self.atoms, dimension=1)
+
+            # gr_labels  = extract_from_list(b, self.labels, dimension=1)
+            # gr_coord   = extract_from_list(b, self.coord, dimension=1)
+            # gr_radii   = extract_from_list(b, self.radii, dimension=1)
+            # gr_atoms   = extract_from_list(b, self.atoms, dimension=1)
+            print(gr_indices, b, conn_idx, self.labels, len(self.labels), gr_labels, gr_radii)
+            print(len(self.atoms))
+            print(gr_atoms)
+            newgroup = group(gr_labels, gr_coord, parent_indices=gr_indices, radii=gr_radii, parent=self)
             newgroup.set_atoms(atomlist=gr_atoms, overwrite_parent=False)
             newgroup.get_closest_metal()
             self.groups.append(newgroup)
-
+        print(self.groups)
+
     #######################################################
     def get_hapticity(self, debug: int=0):
         if not hasattr(self,"groups"): self.split_ligand(debug=debug)
@@ -540,13 +551,23 @@ def __init__(self, labels: list, coord: list, parent_indices: list=None, radii:
         specie.__init__(self, labels, coord, parent_indices, radii, parent)
 
     #######################################################
+    # def __repr__(self):
+    #     to_print = ""
+    #     to_print += f'---------------------------------------------------\n'
+    #     to_print +=  '   >>> Cell2mol GROUP Object >>>                   \n'
+    #     to_print += f'---------------------------------------------------\n'
+    #     specie.__repr__(self, indirect=True)
     def __repr__(self):
-        to_print = ""
-        to_print += f'---------------------------------------------------\n'
-        to_print +=  '   >>> Cell2mol GROUP Object >>>                   \n'
+        to_print  = f'---------------------------------------------------\n'
+        to_print +=  '   >>> Cell2mol GROUP Object >>>                    \n'
         to_print += f'---------------------------------------------------\n'
-        specie.__repr__(self, indirect=True)
-
+        to_print += f' Version               = {self.version}\n'
+        to_print += f' Type                  = {self.type}\n'
+        if hasattr(self,'subtype'): to_print += f' Sub-Type              = {self.subtype}\n'
+        to_print += f' Label                 = {self.labels}\n'
+        to_print += f' Indices in Parent       = {self.parent_indices}\n'
+        to_print += '----------------------------------------------------\n'
+        return to_print
     #######################################################
     def get_closest_metal(self, debug: int=0):
         apos = compute_centroid(np.array(self.coord))
@@ -585,13 +606,13 @@ def get_hapticity(self, debug: int=0):
         return self.haptic_type 
 
     #######################################################
-    def coordination_correction(self, debug=1) :
+    def coordination_correction(self, debug=2) :
         if self.is_haptic:                     self = coordination_correction_for_haptic(self, debug=debug)
         if self.is_haptic == False :           self = coordination_correction_for_nonhaptic(self, debug=debug)
         return self 
 
     #######################################################
-    def check_denticity(self, debug: int=0):
+    def check_denticity(self, debug: int=2):
         from cell2mol.connectivity import add_atom
         if not hasattr(self,"is_haptic"): self.get_hapticity()
         if not hasattr(self,"atoms"):     self.set_atoms()
@@ -720,14 +741,14 @@ def get_connected_metals(self, metalist: list, debug: int=0):
         return self.metals
 
     #######################################################
-    def get_closest_metal(self, metalist: list, debug: int=0):
+    def get_closest_metal(self, debug: int=0):
         ## Here, the list of metal atoms must be provided
         apos = self.coord
         dist = []
-        for met in metalist:
+        for met in self.parent.parent.metals:
             bpos = np.array(met.coord)
             dist.append(np.linalg.norm(apos - bpos))
-        self.closest_metal = metalist[np.argmin(dist)]
+        self.closest_metal = self.parent.parent.metals[np.argmin(dist)]
         return self.closest_metal
 
     #######################################################
@@ -753,6 +774,7 @@ def __repr__(self):
 
     #######################################################
     def reset_mconnec(self, idx: int):
+        print(f"ATOM.RESET_MCONN: resetting mconnec for atom {self.label=} with {self.parent_index=} {self.parent=} index {idx=}")
         self.mconnec = 0                                                                        # Corrects data of atom object in self
         self.parent.atoms[self.parent_indices[idx]].mconnec = 0                                 # Corrects data of atom object in ligand class
         self.parent.madjnum[self.parent_indices[idx]] = 0                                       # Corrects data in metal_adjacency number of the ligand class
@@ -970,7 +992,7 @@ def get_reference_molecules(self, ref_labels: list, ref_fracs: list, cov_factor:
                     ref.get_hapticity(debug=debug)                          ### Former "get_hapticity(ref)" function
                     # ref.get_coordination_geometry(debug=debug)                ### Former "get_coordination_Geometry(ref)" function 
                     for lig in ref.ligands:
-                        lig.get_denticity(debug=0)
+                        lig.get_denticity(debug=2)
 
         if isgood: self.has_isolated_H = False
         else:      self.has_isolated_H = True

cell2mol/coordination_sphere.py

@@ -318,6 +318,7 @@ def coordination_correction_for_nonhaptic (group, debug=1) -> list:
     ## Second Correction
     for idx, atom in enumerate(group.atoms):
         metal = atom.get_closest_metal()
+        dist = get_dist(atom.coord, metal.coord)
         thres = get_thres_from_two_atoms(metal.label, atom.label, debug=debug)
         if debug >= 1 : print(f"\tAtom {atom.label} connected to {metal.label} distance {get_dist(atom.coord, metal.coord)} with threshold {thres}")
 
@@ -369,27 +370,31 @@ def coordination_correction_for_nonhaptic (group, debug=1) -> list:
     return group 
 
 #######################################################    
-def coordination_correction_for_haptic (group, debug=1) -> list:
+def coordination_correction_for_haptic (group, debug=2) -> list:
 
     ratio_list = []
     for idx, atom in enumerate(group.atoms):
         metal = atom.get_closest_metal()
         dist = get_dist(atom.coord, metal.coord)
         thres = get_thres_from_two_atoms(metal.label, atom.label, debug=debug)
         ratio_list.append(round(dist/thres,3))
-        if debug >= 1 : print(f"\tAtom {idx} :", atom.label, f"\tMetal :", metal.label, "\tdistance :", round(dist, 3), "\tthres :", thres)
+        # if debug >= 1 : 
+        print(f"\tAtom {idx} :", atom.label, f"\tMetal :", metal.label, "\tdistance :", round(dist, 3), "\tthres :", thres)
 
     std_dev = round(np.std(ratio_list), 3)
-    if debug >= 1 : print(f"{ratio_list=} {std_dev=}")
-
+    # if debug >= 1 : 
+    print(f"{ratio_list=} {std_dev=}")
+    # print()
     count = 0
     for idx, (atom, ratio) in enumerate(zip(group.atoms, ratio_list)) :
         if atom.label == "H" : 
             if debug >=1 : print("\t!!! Wrong metal-coordination assignment for Atom", idx, atom.label , get_dist(atom.coord, metal.coord), "due to H")
+            print(atom.label)
             atom.reset_mconnec(idx)
             count += 1          
         elif std_dev > 0.05 and ratio > 0.9 :
             if debug >=1 : print("\t!!! Wrong metal-coordination assignment for Atom", idx, atom.label , get_dist(atom.coord, metal.coord), "due to the long distance")
+            print(atom.label)
             atom.reset_mconnec(idx) 
             count += 1      
         else :

cell2mol/other.py

@@ -1,6 +1,6 @@
 ## File for small functions
 import numpy as np
-
+import sys
 ################################
 def extract_from_list(entrylist: list, old_array: list, dimension: int=2, debug: int=0) -> list:
     #if debug >= 0: print("EXTRACT_FROM_LIST. received:", len(entrylist), np.max(entrylist)+1, len(old_array))