Merge pull request #205 from hjkgrp/first_shell

added new function to obtain first coordination shell while accountin…

molSimplify/Classes/mol3D.py

@@ -30,7 +30,6 @@
 try:
     import PyQt5  # noqa: F401
     from molSimplify.Classes.miniGUI import miniGUI
-
     # PyQt5 flag
     qtflag = True
 except ImportError:
@@ -1185,6 +1184,65 @@ def get_smilesOBmol_charge(self):
         charge = charge - nh_obmol + nh
         return charge
 
+    def get_first_shell(self, check_hapticity=True):
+        '''
+        Get the first coordination shell of a mol3D object with a single transition metal (read from CSD mol2 file)
+        if check_hapticity is True updates the first shell of multiheptate ligand to be hydrogen set at the geometric mean
+
+        Parameters
+        ----------
+            check_hapticity: boolean
+                whether to update multiheptate ligands to their geometric centroid
+        Returns
+        ----------
+            mol 3D object: first coordination shell with metal (can change based on check_hapticity)
+            list: list of hapticity
+        '''
+        from molSimplify.Informatics.graph_analyze import obtain_truncation_metal
+        import networkx as nx
+        mol_fcs = obtain_truncation_metal(self, hops=1)
+        M_coord = mol_fcs.getAtomCoords(mol_fcs.findMetal()[0])
+        M_sym = mol_fcs.getAtom(mol_fcs.findMetal()[0]).symbol()
+        G = nx.from_numpy_array(mol_fcs.graph)
+        G.remove_node(mol_fcs.findMetal()[0])
+        coord_list = [c for c in sorted(nx.connected_components(G), key=len, reverse=True)]
+        hapticity_list = [len(c) for c in sorted(nx.connected_components(G), key=len, reverse=True)]
+        new_coords_mol = []
+        new_coords_sym = []
+        if not len(coord_list) == G.number_of_nodes():
+            for i in range(len(coord_list)):
+                if len(coord_list[i]) == 1:
+                    coord_index = list(coord_list[i])[0]
+                    coord = mol_fcs.getAtomCoords(coord_index)
+                    sym = mol_fcs.getAtom(coord_index).symbol()
+                    new_coords_mol.append(coord)
+                    new_coords_sym.append(sym)
+                else:
+                    get_centroid = []
+                    for j in coord_list[i]:
+                        get_centroid.append(mol_fcs.getAtomCoords(j))
+                    coordinating = np.array(get_centroid)
+                    coord = np.mean(coordinating, axis=0)
+                    new_coords_mol.append(coord.tolist())
+                    new_coords_sym.append('H')
+            new_mol = mol3D()
+            new_mol.bo_dict = {}
+            new_mol.addAtom(atom3D(M_sym, M_coord))
+            for i in range(len(new_coords_mol)):
+                new_mol.addAtom(atom3D(new_coords_sym[i], new_coords_mol[i]))
+            new_mol.graph = np.zeros([new_mol.natoms, new_mol.natoms])
+            for i in range(new_mol.natoms):
+                if i != new_mol.findMetal()[0]:
+                    new_mol.add_bond(new_mol.findMetal()[0], i, 1)
+        else:
+            new_mol = mol3D()
+            new_mol.copymol3D(mol_fcs)
+
+        if check_hapticity:
+            return new_mol, hapticity_list
+        else:
+            return mol_fcs, hapticity_list
+
     def get_octetrule_charge(self, debug=False):
         '''
         Get the octet-rule charge provided a mol3D object with bo_graph (read from CSD mol2 file)
@@ -5385,7 +5443,7 @@ def is_edge_compound(self, transition_metals_only: bool = True) -> Tuple[int, Li
             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=None,
+    def get_geometry_type_old(self, dict_check=False, angle_ref=False, num_coord=None,
                           flag_catoms=False, catoms_arr=None, debug=False,
                           skip=False, transition_metals_only=False, num_recursions=[0, 0]):
         """
@@ -5463,7 +5521,7 @@ def get_geometry_type(self, dict_check=False, angle_ref=False, num_coord=None,
                 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])
+            return mol_copy.get_geometry_type_old(num_recursions=[num_sandwich_lig, num_edge_lig])
 
         if num_edge_lig:
             mol_copy = mol3D()
@@ -5482,7 +5540,7 @@ def get_geometry_type(self, dict_check=False, angle_ref=False, num_coord=None,
                 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])
+            return mol_copy.get_geometry_type_old(num_recursions=[num_sandwich_lig, num_edge_lig])
 
         if num_coord not in all_geometries:
             geometry = "unknown"
@@ -5517,9 +5575,111 @@ def get_geometry_type(self, dict_check=False, angle_ref=False, num_coord=None,
             "angle_devi": angle_devi,
             "summary": summary,
             "num_sandwich_lig": num_recursions[0],
+            "info_sandwich_lig": info_sandwich_lig,
             "aromatic": aromatic,
             "allconnect": allconnect,
-            "num_edge_lig": num_recursions[1]
+            "num_edge_lig": num_recursions[1],
+            "info_edge_lig": info_edge_lig,
+        }
+        return results
+
+    def get_geometry_type(self, dict_check=False, angle_ref=False,
+                          flag_catoms=False, catoms_arr=None, debug=False,
+                          skip=False, transition_metals_only=False):
+        """
+        Get the type of the geometry (linear (2), trigonal planar(3), tetrahedral(4), square planar(4),
+        trigonal bipyramidal(5), square pyramidal(5, one-empty-site),
+        octahedral(6), pentagonal bipyramidal(7))
+
+        uses hapticity truncated first coordination shell.
+        Does not require the input of num_coord.
+
+        Parameters
+        ----------
+            dict_check : dict, optional
+                The cutoffs of each geo_check metrics we have. Default is False
+            angle_ref : bool, optional
+                Reference list of list for the expected angles (A-metal-B) of each connection atom.
+            num_coord : int, optional
+                Expected coordination number.
+            flag_catoms : bool, optional
+                Whether or not to return the catoms arr. Default as False.
+            catoms_arr : Nonetype, optional
+                Uses the catoms of the mol3D by default. User and overwrite this connection atom array by explicit input.
+                Default is Nonetype.
+            debug : bool, optional
+                Flag for extra printout. Default is False.
+            skip : list, optional
+                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.
+
+        Returns
+        -------
+            results : dictionary
+                Measurement of deviations from arrays.
+
+        """
+
+        first_shell, hapt = self.get_first_shell()
+        num_coord = first_shell.natoms - 1
+        all_geometries = globalvars().get_all_geometries()
+        all_angle_refs = globalvars().get_all_angle_refs()
+        summary = {}
+
+        if len(first_shell.graph):  # Find num_coord based on metal_cn if graph is assigned
+            if len(first_shell.findMetal()) > 1:
+                raise ValueError('Multimetal complexes are not yet handled.')
+            elif len(first_shell.findMetal(transition_metals_only=transition_metals_only)) == 1:
+                num_coord = len(first_shell.getBondedAtomsSmart(first_shell.findMetal(transition_metals_only=transition_metals_only)[0]))
+            else:
+                raise ValueError('No metal centers exist in this complex.')
+
+        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.")
+
+        if num_coord not in [2, 3, 4, 5, 6, 7]:
+            results = {
+                "geometry": "unknown",
+                "angle_devi": False,
+                "summary": {},
+                "hapticity": hapt,
+            }
+            return results
+        elif num_coord == 2:
+            if first_shell.findMetal()[0] == 2:
+                angle = first_shell.getAngle(0, 2, 1)
+            elif first_shell.findMetal()[0] == 1:
+                angle = first_shell.getAngle(0, 1, 2)
+            else:
+                angle = first_shell.getAngle(1, 0, 2)
+            results = {
+                "geometry": "linear",
+                "angle_devi": 180 - angle,
+                "summary": {},
+                "hapticity": hapt,
+            }
+            return results
+
+        possible_geometries = all_geometries[num_coord]
+        for geotype in possible_geometries:
+            dict_catoms_shape, catoms_assigned = first_shell.oct_comp(
+                angle_ref=all_angle_refs[geotype], catoms_arr=None, debug=debug)
+            if debug:
+                print("Geocheck assigned catoms: ", catoms_assigned,
+                      [first_shell.getAtom(ind).symbol() for ind in catoms_assigned])
+            summary.update({geotype: dict_catoms_shape})
+
+        angle_devi, geometry = 10000, None
+        for geotype in summary:
+            if summary[geotype]["oct_angle_devi_max"] < angle_devi:
+                angle_devi = summary[geotype]["oct_angle_devi_max"]
+                geometry = geotype
+        results = {
+            "geometry": geometry,
+            "angle_devi": angle_devi,
+            "summary": summary,
+            "hapticity": hapt,
         }
         return results
 

tests/test_mol3D.py

@@ -106,34 +106,32 @@ def test_mutating_atoms():
     assert mol.findMetal() == []
 
 
-@pytest.mark.parametrize('name, coordination_number, geometry_str', [
-    ('linear', 2, 'linear'),
-    ('trigonal_planar', 3, 'trigonal planar'),
-    ('t_shape', 3, 'T shape'),
-    ('trigonal_pyramidal', 3, 'trigonal pyramidal'),
-    ('tetrahedral', 4, 'tetrahedral'),
-    ('square_planar', 4, 'square planar'),
-    ('seesaw', 4, 'seesaw'),
-    ('trigonal_bipyramidal', 5, 'trigonal bipyramidal'),
-    ('square_pyramidal', 5, 'square pyramidal'),
-    # ('pentagonal_planar', 5, 'pentagonal planar'),
-    ('octahedral', 6, 'octahedral'),
-    # ('pentagonal_pyramidal', 6, 'pentagonal pyramidal'),
-    ('trigonal_prismatic', 6, 'trigonal prismatic'),
-    # ('pentagonal_bipyramidal', 7, 'pentagonal bipyramidal')
-    # ('square_antiprismatic', 8, 'square antiprismatic'),
-    # ('tricapped_trigonal_prismatic', 9, 'tricapped trigonal prismatic'),
+@pytest.mark.parametrize('name, geometry_str', [
+    ('linear', 'linear'),
+    ('trigonal_planar', 'trigonal planar'),
+    ('t_shape', 'T shape'),
+    ('trigonal_pyramidal', 'trigonal pyramidal'),
+    ('tetrahedral', 'tetrahedral'),
+    ('square_planar', 'square planar'),
+    ('seesaw', 'seesaw'),
+    ('trigonal_bipyramidal', 'trigonal bipyramidal'),
+    ('square_pyramidal', 'square pyramidal'),
+    # ('pentagonal_planar', 'pentagonal planar'),
+    ('octahedral', 'octahedral'),
+    # ('pentagonal_pyramidal', 'pentagonal pyramidal'),
+    ('trigonal_prismatic', 'trigonal prismatic'),
+    # ('pentagonal_bipyramidal', 'pentagonal bipyramidal')
+    # ('square_antiprismatic', 'square antiprismatic'),
+    # ('tricapped_trigonal_prismatic', 'tricapped trigonal prismatic'),
 ])
-def test_get_geometry_type(resource_path_root, name, coordination_number, geometry_str):
+def test_get_geometry_type(resource_path_root, name, geometry_str):
     xyz_file = resource_path_root / "inputs" / "geometry_type" / f"{name}.xyz"
     mol = mol3D()
     mol.readfromxyz(xyz_file)
 
-    geo_report = mol.get_geometry_type(num_coord=coordination_number, debug=True)
+    geo_report = mol.get_geometry_type(debug=True)
 
     assert geo_report['geometry'] == geometry_str
-    assert geo_report['allconnect'] is False
-    assert geo_report['aromatic'] is False
 
 
 def test_get_geometry_type_catoms_arr(resource_path_root):
@@ -142,34 +140,38 @@ def test_get_geometry_type_catoms_arr(resource_path_root):
     mol.readfromxyz(xyz_file)
 
     with pytest.raises(ValueError):
-        mol.get_geometry_type(num_coord=6, catoms_arr=[1], debug=True)
+        mol.get_geometry_type(catoms_arr=[1], debug=True)
 
-    geo_report = mol.get_geometry_type(num_coord=6, catoms_arr=[1, 4, 7, 10, 13, 16], debug=True)
+    geo_report = mol.get_geometry_type(catoms_arr=[1, 4, 7, 10, 13, 16], debug=True)
 
     assert geo_report['geometry'] == 'octahedral'
-    assert geo_report['allconnect'] is False
-    assert geo_report['aromatic'] is False
 
 
 @pytest.mark.parametrize(
-    'name, geometry_str, num_sandwich',
+    'name, geometry_str, hapticity',
     [
-        ('BOWROX_comp_0.mol2', 'tetrahedral', 1),
-        ('BOXTEQ_comp_0.mol2', 'tetrahedral', 1),
-        ('BOXTIU_comp_0.mol2', 'tetrahedral', 1),
-        # ('BOZHOQ_comp_2.mol2', 'linear', 2),
-        # ('BOZHUW_comp_2.mol2', 'linear', 2),
+        ('BOWROX_comp_0.mol2', 'tetrahedral', [5, 1, 1, 1]),
+        ('BOXTEQ_comp_0.mol2', 'tetrahedral', [6, 1, 1, 1]),
+        ('BOXTIU_comp_0.mol2', 'tetrahedral', [6, 1, 1, 1]),
+        ('BOZHOQ_comp_2.mol2', 'linear', [5, 5]),
+        ('BOZHUW_comp_2.mol2', 'linear', [5, 5]),
+        ('BUFLUM_comp_0.mol2', 'T shape', [2, 1, 1]),
+        ('BUHMID_comp_0.mol2', 'trigonal planar', [3, 1, 1]),
+        ('COYXUM_comp_0.mol2', 'tetrahedral', [5, 1, 1, 1]),
+        ('COYYEX_comp_0.mol2', 'trigonal planar', [5, 1, 1]),
+        ('COYYIB_comp_0.mol2', 'tetrahedral', [5, 1, 1, 1]),
     ]
 )
-def test_get_geometry_type_sandwich(resource_path_root, name, geometry_str, num_sandwich):
-    input_file = resource_path_root / "inputs" / "sandwich_compounds" / name
+def test_get_geometry_type_hapticity(resource_path_root, name, geometry_str, hapticity):
+    input_file = resource_path_root / "inputs" / "hapticity_compounds" / name
     mol = mol3D()
     mol.readfrommol2(input_file)
 
     geo_report = mol.get_geometry_type(debug=True)
 
+    print(geo_report)
     assert geo_report["geometry"] == geometry_str
-    assert geo_report["num_sandwich_lig"] == num_sandwich
+    assert geo_report["hapticity"] == hapticity
 
 
 @pytest.mark.parametrize(
@@ -183,7 +185,7 @@ def test_get_geometry_type_sandwich(resource_path_root, name, geometry_str, num_
     ]
 )
 def test_is_sandwich_compound(resource_path_root, name, con_atoms):
-    input_file = resource_path_root / "inputs" / "sandwich_compounds" / name
+    input_file = resource_path_root / "inputs" / "hapticity_compounds" / name
     mol = mol3D()
     mol.readfrommol2(input_file)
 
@@ -204,13 +206,10 @@ def test_is_sandwich_compound(resource_path_root, name, con_atoms):
     [
         ("BUFLUM_comp_0.mol2", [{2, 4}]),
         ("BUHMID_comp_0.mol2", [{3, 4, 5}]),
-        ("COYXUM_comp_0.mol2", [{0, 1, 2, 3, 4}]),
-        ("COYYEX_comp_0.mol2", [{0, 1, 2, 3, 4}]),
-        ("COYYIB_comp_0.mol2", [{1, 2, 3, 6, 7}]),
     ]
 )
 def test_is_edge_compound(resource_path_root, name, con_atoms):
-    input_file = resource_path_root / "inputs" / "edge_compounds" / name
+    input_file = resource_path_root / "inputs" / "hapticity_compounds" / name
     mol = mol3D()
     mol.readfrommol2(input_file)