Updated get_geometry to properly classify edge / sandwich ligands and…

… coordination numbers > 7

molSimplify/Classes/globalvars.py

@@ -402,27 +402,46 @@
                                  'devi_linear_avrg': 20, 'devi_linear_max': 28}
                        }
 
+dict_eightcoord_check = {"mono": {'num_coord_metal': 8,
+                              'rmsd_max': 0.4, 'atom_dist_max': 0.6,
+                              'oct_angle_devi_max': 16, 'max_del_sig_angle': 27,
+                              'dist_del_eq': 0.45, 'dist_del_all': 1.25,
+                              'devi_linear_avrg': 35, 'devi_linear_max': 40},
+                         "multi": {'num_coord_metal': 8,
+                               'rmsd_max': 0.4, 'atom_dist_max': 0.6,
+                               'oct_angle_devi_max': 20, 'max_del_sig_angle': 35,
+                               'dist_del_eq': 0.45, 'dist_del_all': 1.25,
+                               'devi_linear_avrg': 35, 'devi_linear_max': 40}
+                         }
+
 dict_staus = {'good': 1, 'bad': 0}
 
 oct_angle_ref = [[90, 90, 90, 90, 180] for x in range(6)]
 tetra_angle_ref = [[109.47, 109.47, 109.47] for x in range(4)]
 oneempty_angle_ref = [[90, 90, 90, 90], [180, 90, 90, 90], [180, 90, 90, 90],
                       [180, 90, 90, 90], [180, 90, 90, 90]]
+eightcoord_angle_ref = [[73.7, 73.7, 76.8, 76.8, 116.1, 142.1, 142.1] for x in range(8)]
 geo_check_dictionary = {"dict_oct_check_loose": dict_oct_check_loose,
                         "dict_oct_check_st": dict_oct_check_st,
                         "dict_oneempty_check_st": dict_oneempty_check_st,
                         "dict_oneempty_check_loose": dict_oneempty_check_loose,
                         "dict_staus": dict_staus,
                         "oct_angle_ref": oct_angle_ref,
-                        "oneempty_angle_ref": oneempty_angle_ref}
+                        "oneempty_angle_ref": oneempty_angle_ref,
+                        "eightcoord_angle_ref": eightcoord_angle_ref}
 all_geometries = {
+    2: ["linear", "bent"],
     3: ["trigonal planar", "T shape", "trigonal pyramidal"],
     4: ["tetrahedral", "square planar", "seesaw"],
     5: ["trigonal bipyramidal", "square pyramidal", "pentagonal planar"],
     6: ["octahedral", "pentagonal pyramidal", "trigonal prismatic"],
     7: ["pentagonal bipyramidal"],
+    8: ["square antiprismatic"],
+    9: ["tricapped trigonal prismatic"]
 }
 all_angle_refs = {
+    "linear": [[180] for x in range(2)],
+    "bent": [[120] for x in range(2)],
     "trigonal planar": [[120, 120] for x in range(3)],
     "T shape": [[90, 90], [90, 180], [90, 180]],
     "tetrahedral": [[109.47, 109.47, 109.47] for x in range(4)],
@@ -435,10 +454,14 @@
     "octahedral": [[90, 90, 90, 90, 180] for x in range(6)],
     "pentagonal pyramidal": [[90, 90, 90, 90, 90]] + [[36, 36, 72, 72, 90] for x in range(5)],
     "trigonal prismatic": [[75.3, 86.5, 86.5, 133.3, 133.3] for x in range(6)],
-    "pentagonal bipyramidal": [[90, 90, 90, 90, 90, 180] for x in range(2)] + [[72, 72, 144, 144, 90, 90] for x in range(5)]
+    "pentagonal bipyramidal": [[90, 90, 90, 90, 90, 180] for x in range(2)] +
+                              [[72, 72, 144, 144, 90, 90] for x in range(5)],
+    "square antiprismatic": [[73.7, 73.7, 76.8, 76.8, 116.1, 142.1, 142.1] for x in range(8)],
+    "tricapped trigonal prismatic": [[70.5, 70.5, 120, 67.7, 67.7, 135.5, 135.5, 120] for x in range(4)] +
+                                    [[76.3, 79, 138.2, 70.5, 67.7, 76.3, 138.2, 135.5] for x in range(4)] +
+                                    [[70.5, 70.5, 67.7, 67.7, 120, 120, 135.5, 135.5]]
 }
 
-
 # Module for running bash commands
 #  @param cmd String containing command to be run
 #  @return bash output string

molSimplify/Classes/mol3D.py

@@ -5299,12 +5299,14 @@ def is_sandwich_compound(self, transition_metals_only=True):
                         _sl.append(set(lig))
                         break
         num_sandwich_lig = len(sandwich_ligands)
-        info_sandwich_lig = [{"natoms_connected": len(
-            x[0]), "natoms_ring": x[1], "aromatic": x[2]} for x in sandwich_ligands]
+        info_sandwich_lig = [
+                {"natoms_connected": len(x[0]), "natoms_ring": x[1], "aromatic": x[2]} for x in sandwich_ligands]
+        sandwich_lig_atoms = [
+                {"atom_idxs": x[0]} for x in sandwich_ligands]
         aromatic = any([x["aromatic"] for x in info_sandwich_lig])
         allconnect = any([x["natoms_connected"] == x["natoms_ring"]
                           for x in info_sandwich_lig])
-        return num_sandwich_lig, info_sandwich_lig, aromatic, allconnect
+        return num_sandwich_lig, info_sandwich_lig, aromatic, allconnect, list(sandwich_lig_atoms)
 
     def is_edge_compound(self, transition_metals_only=True):
         """
@@ -5323,7 +5325,7 @@ def is_edge_compound(self, transition_metals_only=True):
         # two connected non-metal atoms both connected to the metal.
 
         from molSimplify.Informatics.graph_analyze import obtain_truncation_metal
-        num_sandwich_lig, info_sandwich_lig, aromatic, allconnect = self.is_sandwich_compound(transition_metals_only=transition_metals_only)
+        num_sandwich_lig, info_sandwich_lig, aromatic, allconnect, sandwich_lig_atoms = self.is_sandwich_compound(transition_metals_only=transition_metals_only)
         if not num_sandwich_lig or (num_sandwich_lig and not allconnect):
             mol_fcs = obtain_truncation_metal(self, hops=1)
             metal_ind = mol_fcs.findMetal(transition_metals_only=transition_metals_only)[0]
@@ -5336,17 +5338,19 @@ def is_edge_compound(self, transition_metals_only=True):
                 if len(lig) >= 2 and not set(lig) in _el:
                     edge_ligands.append([set(lig)])
                     _el.append(set(lig))
-                    break
+                    # break
             num_edge_lig = len(edge_ligands)
             info_edge_lig = [
                 {"natoms_connected": len(x[0])} for x in edge_ligands]
+            edge_lig_atoms = [
+                {"atom_idxs": x[0]} for x in edge_ligands]
         else:
-            num_edge_lig, info_edge_lig = 0, list()
-        return num_edge_lig, info_edge_lig
+            num_edge_lig, info_edge_lig, edge_lig_atoms = 0, list(), list()
+        return num_edge_lig, info_edge_lig, edge_lig_atoms
 
     def get_geometry_type(self, dict_check=False, angle_ref=False, num_coord=False,
                           flag_catoms=False, catoms_arr=None, debug=False,
-                          skip=False, transition_metals_only=False):
+                          skip=False, transition_metals_only=False, num_recursions=[0, 0]):
         """
         Get the type of the geometry (trigonal planar(3), tetrahedral(4), square planar(4),
         trigonal bipyramidal(5), square pyramidal(5, one-empty-site),
@@ -5371,13 +5375,16 @@ def get_geometry_type(self, dict_check=False, angle_ref=False, num_coord=False,
                 Geometry checks to skip. Default is False.
             transition_metals_only : bool, optional
                 Flag for considering more than just transition metals as metals. Default is False.
+            num_recursions : list, optional
+                counter to track number of ligands classified as 'sandwich' and 'edge' in original structure
 
         Returns
         -------
             results : dictionary
                 Measurement of deviations from arrays.
 
         """
+        # from molSimplify.Classes.ligand import ligand_breakdown
 
         all_geometries = globalvars().get_all_geometries()
         all_angle_refs = globalvars().get_all_angle_refs()
@@ -5388,7 +5395,6 @@ def get_geometry_type(self, dict_check=False, angle_ref=False, num_coord=False,
                 raise ValueError('Multimetal complexes are not yet handled.')
             elif len(self.findMetal(transition_metals_only=transition_metals_only)) == 1:
                 num_coord = len(self.getBondedAtomsSmart(self.findMetal(transition_metals_only=transition_metals_only)[0]))
-                # print("coord number:", num_coord)
             else:
                 raise ValueError('No metal centers exist in this complex.')
 
@@ -5400,26 +5406,57 @@ def get_geometry_type(self, dict_check=False, angle_ref=False, num_coord=False,
         if catoms_arr is not None and len(catoms_arr) != num_coord:
             raise ValueError("num_coord and the length of catoms_arr do not match.")
 
-        num_sandwich_lig, info_sandwich_lig, aromatic, allconnect = self.is_sandwich_compound(transition_metals_only=transition_metals_only)
-        num_edge_lig, info_edge_lig = self.is_edge_compound(transition_metals_only=transition_metals_only)
+        num_sandwich_lig, info_sandwich_lig, aromatic, allconnect, sandwich_lig_atoms = self.is_sandwich_compound(transition_metals_only=transition_metals_only)
+        num_edge_lig, info_edge_lig, edge_lig_atoms = self.is_edge_compound(transition_metals_only=transition_metals_only)
 
-        if num_coord not in [3, 4, 5, 6, 7]:
-            if num_sandwich_lig:
-                geometry = "sandwich"
-            elif num_edge_lig:
-                geometry = "edge"
-            else:
-                geometry = "unknown"
+        if num_sandwich_lig:
+            mol_copy = mol3D()
+            mol_copy.copymol3D(mol0=self)
+            catoms = mol_copy.getBondedAtoms(idx=self.findMetal()[0])
+            centroid_coords = []
+            sandwich_lig_catom_idxs = []
+            for idx in range(num_sandwich_lig):
+                sandwich_lig_catoms = np.array(list(sandwich_lig_atoms[idx]['atom_idxs']))-1
+                sandwich_lig_catom_idxs.extend([catoms[i] for i in sandwich_lig_catoms])
+                atom_coords = np.array([mol_copy.getAtomCoords(idx=atom_idx) for atom_idx in [catoms[i] for i in sandwich_lig_catoms]])
+                centroid_coords.append([np.mean(atom_coords[:, 0]), np.mean(atom_coords[:, 1]), np.mean(atom_coords[:, 2])])
+            mol_copy.deleteatoms(sandwich_lig_catom_idxs)
+            for idx in range(num_sandwich_lig):
+                atom = atom3D()
+                atom.setcoords(xyz=centroid_coords[idx])
+                mol_copy.addAtom(atom)
+                mol_copy.add_bond(idx1=mol_copy.findMetal()[0], idx2=mol_copy.natoms-1, bond_type=1)
+            return mol_copy.get_geometry_type(num_recursions=[num_sandwich_lig, num_edge_lig])
+
+        if num_edge_lig:
+            mol_copy = mol3D()
+            mol_copy.copymol3D(mol0=self)
+            catoms = mol_copy.getBondedAtoms(idx=self.findMetal()[0])
+            centroid_coords = []
+            edge_lig_catom_idxs = []
+            for idx in range(num_edge_lig):
+                edge_lig_catoms = np.array(list(edge_lig_atoms[idx]['atom_idxs']))-1
+                edge_lig_catom_idxs.extend([catoms[i] for i in edge_lig_catoms])
+                atom_coords = np.array([mol_copy.getAtomCoords(idx=atom_idx) for atom_idx in [catoms[i] for i in edge_lig_catoms]])
+                centroid_coords.append([np.mean(atom_coords[:, 0]), np.mean(atom_coords[:, 1]), np.mean(atom_coords[:, 2])])
+            mol_copy.deleteatoms(edge_lig_catom_idxs)
+            for idx in range(num_edge_lig):
+                atom = atom3D()
+                atom.setcoords(xyz=centroid_coords[idx])
+                mol_copy.addAtom(atom)
+                mol_copy.add_bond(idx1=mol_copy.findMetal()[0], idx2=mol_copy.natoms-1, bond_type=1)
+            return mol_copy.get_geometry_type(num_recursions=[num_sandwich_lig, num_edge_lig])
+
+        if num_coord not in [2, 3, 4, 5, 6, 7, 8, 9]:
+            geometry = "unknown"
             results = {
                 "geometry": geometry,
                 "angle_devi": False,
                 "summary": {},
-                "num_sandwich_lig": num_sandwich_lig,
-                "info_sandwich_lig": info_sandwich_lig,
+                "num_sandwich_lig": num_recursions[0],
                 "aromatic": aromatic,
                 "allconnect": allconnect,
-                "num_edge_lig": num_edge_lig,
-                "info_edge_lig": info_edge_lig,
+                "num_edge_lig": num_recursions[1]
             }
             return results
 
@@ -5438,22 +5475,14 @@ def get_geometry_type(self, dict_check=False, angle_ref=False, num_coord=False,
             if summary[geotype]["oct_angle_devi_max"] < angle_devi:
                 angle_devi = summary[geotype]["oct_angle_devi_max"]
                 geometry = geotype
-        if num_sandwich_lig:
-            geometry = "sandwich"
-            angle_devi = False
-        elif num_edge_lig:
-            geometry = "edge"
-            angle_devi = False
         results = {
             "geometry": geometry,
             "angle_devi": angle_devi,
             "summary": summary,
-            "num_sandwich_lig": num_sandwich_lig,
-            "info_sandwich_lig": info_sandwich_lig,
+            "num_sandwich_lig": num_recursions[0],
             "aromatic": aromatic,
             "allconnect": allconnect,
-            "num_edge_lig": num_edge_lig,
-            "info_edge_lig": info_edge_lig,
+            "num_edge_lig": num_recursions[1]
         }
         return results