diff --git a/molSimplify/Classes/globalvars.py b/molSimplify/Classes/globalvars.py index 54e45098..6ef27f9c 100644 --- a/molSimplify/Classes/globalvars.py +++ b/molSimplify/Classes/globalvars.py @@ -14,6 +14,7 @@ import subprocess from typing import Dict, Tuple from molSimplify.utils.metaclasses import Singleton +import numpy as np # Dictionary containing atomic mass, atomic number, covalent radius, num valence electrons # Data from http://www.webelements.com/ (last accessed May 13th 2015) @@ -463,6 +464,28 @@ [[70.5, 70.5, 67.7, 67.7, 120, 120, 135.5, 135.5]] } +#need each point to be distance 1 from the origin +deg_to_rad = 2*np.pi / 360 +all_polyhedra = { + "linear": np.array([(1, 0, 0), (-1, 0, 0)]), + "bent": np.array([(1, 0, 0), (np.cos(120*deg_to_rad), np.sin(120*deg_to_rad), 0)]), + "trigonal planar": np.array([(1, 0, 0), (np.cos(120*deg_to_rad), np.sin(120*deg_to_rad), 0), (np.cos(240*deg_to_rad), np.sin(240*deg_to_rad), 0)]), + "T shape": np.array([(1, 0, 0), (np.cos(90*deg_to_rad), np.sin(90*deg_to_rad), 0), (np.cos(180*deg_to_rad), np.sin(180*deg_to_rad), 0)]), + "trigonal pyramidal": np.array([(1, -1, -1), (-1, 1, -1), (-1, -1, 1)])/np.sqrt(3), + "tetrahedral": np.array([(1, 1, 1), (1, -1, -1), (-1, 1, -1), (-1, -1, 1)])/np.sqrt(3), + "square planar": np.array([(1, 0, 0), (0, 1, 0), (-1, 0, 0), (0, -1, 0)]), + "seesaw": np.array([(1, 0, 0), (np.cos(120*deg_to_rad), np.sin(120*deg_to_rad), 0), (0, 0, 1), (0, 0, -1)]), + "trigonal bipyramidal": np.array([(0, 0, 1), (0, 0, -1), (1, 0, 0), (np.cos(120*deg_to_rad), np.sin(120*deg_to_rad), 0), (np.cos(240*deg_to_rad), np.sin(240*deg_to_rad), 0)]), + "square pyramidal": np.array([(1, 0, 0), (0, 1, 0), (-1, 0, 0), (0, -1, 0), (0, 0, 1)]), + "pentagonal planar": np.array([(1, 0, 0), (np.cos(72*deg_to_rad), np.sin(72*deg_to_rad), 0), (np.cos(144*deg_to_rad), np.sin(144*deg_to_rad), 0), (np.cos(216*deg_to_rad), np.sin(216*deg_to_rad), 0), (np.cos(288*deg_to_rad), np.sin(288*deg_to_rad), 0)]), + "octahedral": np.array([(1, 0, 0), (0, 1, 0), (-1, 0, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1)]), + "pentagonal pyramidal": np.array([(1, 0, 0), (np.cos(72*deg_to_rad), np.sin(72*deg_to_rad), 0), (np.cos(144*deg_to_rad), np.sin(144*deg_to_rad), 0), (np.cos(216*deg_to_rad), np.sin(216*deg_to_rad), 0), (np.cos(288*deg_to_rad), np.sin(288*deg_to_rad), 0), (0, 0, 1)]), + "trigonal prismatic": np.array([(1, 0, 1), (np.cos(120*deg_to_rad), np.sin(120*deg_to_rad), 1), (np.cos(240*deg_to_rad), np.sin(240*deg_to_rad), 1), (1, 0, -1), (np.cos(120*deg_to_rad), np.sin(120*deg_to_rad), -1), (np.cos(240*deg_to_rad), np.sin(240*deg_to_rad), -1)])/np.sqrt(2), + "pentagonal bipyramidal": np.array([(1, 0, 0), (np.cos(72*deg_to_rad), np.sin(72*deg_to_rad), 0), (np.cos(144*deg_to_rad), np.sin(144*deg_to_rad), 0), (np.cos(216*deg_to_rad), np.sin(216*deg_to_rad), 0), (np.cos(288*deg_to_rad), np.sin(288*deg_to_rad), 0), (0, 0, 1), (0, 0, -1)]), + "square antiprismatic": np.array([(1, 0, 1), (0, 1, 1), (-1, 0, 1), (0, -1, 1), (np.cos(45*deg_to_rad), np.sin(45*deg_to_rad), -1), (np.cos(135*deg_to_rad), np.sin(135*deg_to_rad), -1), (np.cos(225*deg_to_rad), np.sin(225*deg_to_rad), -1), (np.cos(315*deg_to_rad), np.sin(315*deg_to_rad), -1)])/np.sqrt(2), + "tricapped trigonal prismatic": np.array([(1, 0, 1)/np.sqrt(2), (np.cos(120*deg_to_rad), np.sin(120*deg_to_rad), 1)/np.sqrt(2), (np.cos(240*deg_to_rad), np.sin(240*deg_to_rad), 1)/np.sqrt(2), (1, 0, -1)/np.sqrt(2), (np.cos(120*deg_to_rad), np.sin(120*deg_to_rad), -1)/np.sqrt(2), (np.cos(240*deg_to_rad), np.sin(240*deg_to_rad), -1)/np.sqrt(2), (np.cos(60*deg_to_rad), np.sin(60*deg_to_rad), 0), (np.cos(180*deg_to_rad), np.sin(180*deg_to_rad), 0), (np.cos(300*deg_to_rad), np.sin(300*deg_to_rad), 0)]) +} + # Module for running bash commands # @param cmd String containing command to be run # @return bash output string @@ -720,6 +743,16 @@ def get_all_angle_refs(self): """ return all_angle_refs + def get_all_polyhedra(self): + """Get reference polyhedra dict. + + Returns + ------- + all_polyhedra : dict + Reference polyhedra for various geometries. + """ + return all_polyhedra + def add_custom_path(self, path): """Record custom path in ~/.molSimplify file diff --git a/molSimplify/Classes/mol3D.py b/molSimplify/Classes/mol3D.py index 6d629352..9377c43b 100644 --- a/molSimplify/Classes/mol3D.py +++ b/molSimplify/Classes/mol3D.py @@ -16,7 +16,7 @@ from openbabel import openbabel # version 3 style import except ImportError: import openbabel # fallback to version 2 -from typing import List, Optional, Tuple +from typing import List, Optional, Tuple, Dict, Any from scipy.spatial import ConvexHull from molSimplify.utils.decorators import deprecated @@ -25,7 +25,8 @@ from molSimplify.Scripts.geometry import (distance, connectivity_match, vecangle, rotation_params, rotate_around_axis) -from molSimplify.Scripts.rmsd import rigorous_rmsd +from molSimplify.Scripts.rmsd import rigorous_rmsd, kabsch_rmsd, kabsch_rotate +from itertools import permutations try: import PyQt5 # noqa: F401 @@ -4918,7 +4919,7 @@ def Structure_inspection(self, init_mol=None, catoms_arr=None, num_coord=5, dict num_coord=num_coord, debug=debug) return flag_oct, flag_list, dict_oct_info, flag_oct_loose, flag_list_loose - def get_fcs(self, strict_cutoff=False, catom_list=None): + def get_fcs(self, strict_cutoff=False, catom_list=None, max6=True): """ Get first coordination shell of a transition metal complex. @@ -4928,6 +4929,8 @@ def get_fcs(self, strict_cutoff=False, catom_list=None): strict bonding cutoff for fullerene and SACs catom_list : list, optional List of indices of coordinating atoms. + max6 : bool, optional + If True, will return catoms from oct_comp. Returns ------- @@ -4938,7 +4941,7 @@ def get_fcs(self, strict_cutoff=False, catom_list=None): self.get_num_coord_metal(debug=False, strict_cutoff=strict_cutoff, catom_list=catom_list) catoms = self.catoms # print(catoms, [self.getAtom(x).symbol() for x in catoms]) - if len(catoms) > 6: + if max6 and len(catoms) > 6: _, catoms = self.oct_comp(debug=False) fcs = [metalind] + catoms return fcs @@ -5683,6 +5686,171 @@ def get_geometry_type(self, dict_check=False, angle_ref=False, } return results + def get_geometry_type_distance( + self, max_dev=1e6, close_dev=1e-2, + flag_catoms=False, catoms_arr=None, + skip=False, transition_metals_only=False) -> Dict[str, Any]: + """ + Get the type of the geometry (available options in globalvars all_geometries). + + uses hapticity truncated first coordination shell. + Does not require the input of num_coord. + + Parameters + ---------- + max_dev : float, optional + Maximum RMSD allowed between a structure and an ideal geometry before it is classified as unknown. Default is 1e6. + close_dev : float, optional + Maximum difference in RMSD between two classifications allowed before they are compared by maximum single-atom deviation as well. + 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. + 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 + Contains the classified geometry and the RMSD from an ideal structure. + Summary contains a list of the RMSD and the maximum single-atom deviation for all considered geometry types. + + """ + + first_shell, hapt = self.get_first_shell() + num_coord = first_shell.natoms - 1 + all_geometries = globalvars().get_all_geometries() + all_polyhedra = globalvars().get_all_polyhedra() + 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: + # Use oct=False to ensure coordination number based on radius cutoffs only + num_coord = len(first_shell.getBondedAtomsSmart(first_shell.findMetal(transition_metals_only=transition_metals_only)[0], oct=False)) + 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 list(all_geometries.keys()): + # should we indicate somehow that these are unknown due to a different coordination number? + results = { + "geometry": "unknown", + "rmsd": np.NAN, + "summary": {}, + "hapticity": hapt, + "close_rmsds": False + } + return results + + possible_geometries = all_geometries[num_coord] + + # for each same-coordinated geometry, get the minimum RMSD and the maximum single-atom deviation in that pairing + for geotype in possible_geometries: + rmsd_calc, max_dist = self.dev_from_ideal_geometry(all_polyhedra[geotype]) + summary.update({geotype: [rmsd_calc, max_dist]}) + + close_rmsds = False + current_rmsd, geometry = max_dev, "unknown" + for geotype in summary: + # if the RMSD for this structure is the lowest seen so far (within a threshold) + if summary[geotype][0] < (current_rmsd + close_dev): + # if the RMSDs are close, flag this in the summary and classify on second criterion + if np.abs(summary[geotype][0] - current_rmsd) < close_dev: + close_rmsds = True + if summary[geotype][1] < summary[geometry][1]: + # classify based on largest singular deviation + current_rmsd = summary[geotype][0] + geometry = geotype + else: + current_rmsd = summary[geotype][0] + geometry = geotype + + results = { + "geometry": geometry, + "rmsd": current_rmsd, + "summary": summary, + "hapticity": hapt, + "close_rmsds": close_rmsds + } + return results + + def dev_from_ideal_geometry(self, ideal_polyhedron: np.ndarray) -> Tuple[float, float]: + """ + Return the minimum RMSD between a geometry and an ideal polyhedron (with the same average bond distances). + Enumerates all possible indexing of the geometry. As such, only recommended for small systems. + + Parameters + ---------- + ideal_polyhedron: np.array of 3-tuples of coordinates + Reference list of points for an ideal geometry + + Returns + ------- + rmsd: float + Minimum root mean square distance between the fed geometry and the ideal polyhedron + single_dev: float + Maximum distance between any paired points in the fed geometry and the ideal polyhedron. + """ + + metal_idx = self.findMetal() + if len(metal_idx) == 0: + raise ValueError('No metal centers exist in this complex.') + elif len(metal_idx) != 1: + raise ValueError('Multimetal complexes are not yet handled.') + temp_mol = self.get_first_shell()[0] + fcs_indices = temp_mol.get_fcs(max6=False) + # remove metal index from first coordination shell + fcs_indices.remove(temp_mol.findMetal()[0]) + + if len(fcs_indices) != len(ideal_polyhedron): + raise ValueError('The coordination number differs between the two provided structures.') + + # have to redo getting metal_idx with the new mol after running get_first_shell + # want to work with temp_mol since it has the edge and sandwich logic implemented to replace those with centroids + metal_atom = temp_mol.getAtoms()[temp_mol.findMetal()[0]] + fcs_atoms = [temp_mol.getAtoms()[i] for i in fcs_indices] + # construct a np array of the non-metal atoms in the FCS + distances = [] + positions = np.zeros([len(fcs_indices), 3]) + for idx, atom in enumerate(fcs_atoms): + distance = atom.distance(metal_atom) + distances.append(distance) + # shift so the metal is at (0, 0, 0) + positions[idx, :] = np.array(atom.coords()) - np.array(metal_atom.coords()) + + current_min = np.inf + orders = permutations(range(len(ideal_polyhedron))) + max_dist = 0 + + # if desired, make it so the ideal polyhedron has same average bond distance as the mol + # scaled_polyhedron = ideal_polyhedron * np.mean(np.array(distances)) + + # for all possible assignments, find RMSD between ideal and actual structure + ideal_positions = np.zeros([len(fcs_indices), 3]) + for order in orders: + for i in range(len(order)): + # if you wanted to use the same average bond length for all, use the following + # ideal_positions[i, :] = scaled_polyhedron[order[i]] + # if you want to let each ligand scale its length independently, uncomment the following + ideal_positions[i, :] = ideal_polyhedron[order[i]] * distances[i] + rmsd_calc = kabsch_rmsd(ideal_positions, positions) + if rmsd_calc < current_min: + current_min = rmsd_calc + # calculate and store the maximum pairwise distance + rot_ideal = kabsch_rotate(ideal_positions, positions) + diff_matrix = rot_ideal - positions + pairwise_dists = np.sum(diff_matrix**2, axis=1) + max_dist = np.max(pairwise_dists) + + # return minimum RMSD, maximum pairwise distance in that structure + return current_min, max_dist + def get_features(self, lac=True, force_generate=False, eq_sym=False, use_dist=False, NumB=False, Gval=False, size_normalize=False, alleq=False, strict_cutoff=False, catom_list=None, MRdiag_dict={}, depth=3): diff --git a/molSimplify/Data/tbp.dat b/molSimplify/Data/tbp.dat index 32b02f3a..cd5cee6e 100644 --- a/molSimplify/Data/tbp.dat +++ b/molSimplify/Data/tbp.dat @@ -1,6 +1,6 @@ -0.0 0.0 0.0 --1.7865930597 -0.0003546225 -0.0007425232 -1.7865967124 -0.0003546332 -0.0007425118 -2.86499999724e-07 1.382576768 -1.1498771629 -2.92799999801e-07 0.2985354596 1.7728580026 -4.09099999654e-07 -1.6828658296 -0.6318871017 +0.00000 0.00000 0.00000 +2.00000 0.00000 0.00000 +-1.00000 1.73205 0.00000 +-1.00000 -1.73205 0.00000 +0.00000 0.00000 2.00000 +0.00000 0.00000 -2.00000 diff --git a/molSimplify/Data/tpl.dat b/molSimplify/Data/tpl.dat index eaa924a0..74bd8015 100644 --- a/molSimplify/Data/tpl.dat +++ b/molSimplify/Data/tpl.dat @@ -1,4 +1,4 @@ 0.00000 0.00000 0.00000 --1.58487 0.39844 0.52472 -0.77037 -1.41380 0.59470 -0.78682 0.99734 -1.17123 +2.00000 0.00000 0.00000 +-1.00000 1.73205 0.00000 +-1.00000 -1.73205 0.00000 diff --git a/tests/test_mol3D.py b/tests/test_mol3D.py index 507187d4..3e316559 100644 --- a/tests/test_mol3D.py +++ b/tests/test_mol3D.py @@ -2,6 +2,7 @@ import numpy as np from molSimplify.Classes.mol3D import mol3D from molSimplify.Classes.atom3D import atom3D +from molSimplify.Classes.globalvars import globalvars def test_adding_and_deleting_atoms(): @@ -335,3 +336,68 @@ def test_mol3D_from_smiles_benzene(): np.testing.assert_allclose(mol.graph, ref_graph) np.testing.assert_allclose(mol.bo_graph, ref_bo_graph) + + +@pytest.mark.parametrize( + "geo_type, key", + [ + ('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_dev_from_ideal_geometry(resource_path_root, geo_type, key): + mol = mol3D() + mol.readfromxyz(resource_path_root / "inputs" / "geometry_type" / f"{geo_type}.xyz") + + globs = globalvars() + polyhedra = globs.get_all_polyhedra() + rmsd, max_dev = mol.dev_from_ideal_geometry(polyhedra[key]) + + print(polyhedra[key]) + + assert rmsd < 1e-3 + assert max_dev < 1e-3 + + +@pytest.mark.parametrize( + "geo_type, ref", + [ + ('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_geo_geometry_type_distance(resource_path_root, geo_type, ref): + mol = mol3D() + mol.readfromxyz(resource_path_root / "inputs" / "geometry_type" / f"{geo_type}.xyz") + + result = mol.get_geometry_type_distance() + print(result) + assert result['geometry'] == ref diff --git a/tests/testresources/inputs/geometry_type/seesaw.xyz b/tests/testresources/inputs/geometry_type/seesaw.xyz index db118dcb..ed1d9c5d 100644 --- a/tests/testresources/inputs/geometry_type/seesaw.xyz +++ b/tests/testresources/inputs/geometry_type/seesaw.xyz @@ -1,15 +1,15 @@ 13 -07/15/2022 21:19, XYZ structure generated by mol3D Class, molSimplify +04/30/2024 13:23, XYZ structure generated by mol3D Class, molSimplify Fe 0.000000 0.000000 0.000000 -O -2.120000 -0.000421 -0.000881 -H -2.730515 0.026086 0.779322 -H -2.727124 -0.028395 -0.783794 -O 2.120000 -0.000421 -0.000881 -H 2.728823 -0.572514 -0.533990 -H 2.728652 0.572801 0.531233 -O 0.000000 1.629945 -1.355611 -H 0.153088 1.607265 -2.334413 -H -0.153473 2.588291 -1.155199 -O 0.000000 -1.984703 -0.745222 -H 0.126817 -2.825828 -0.236778 -H -0.126096 -2.282875 -1.681767 +O 2.120000 0.000000 0.000000 +H 2.726198 0.027327 -0.782547 +H 2.726184 -0.027327 0.782547 +O -1.060000 1.835974 0.000000 +H -1.410402 2.333646 -0.781116 +H -1.315789 2.388255 0.781116 +O 0.000000 0.000000 2.120000 +H -0.293326 -0.726007 2.726198 +H 0.293326 0.726007 2.726184 +O 0.000000 -0.000000 -2.120000 +H 0.222961 0.750609 -2.726198 +H -0.222961 -0.750609 -2.726184 diff --git a/tests/testresources/inputs/geometry_type/trigonal_bipyramidal.xyz b/tests/testresources/inputs/geometry_type/trigonal_bipyramidal.xyz index 2ffa8c99..55c941ac 100644 --- a/tests/testresources/inputs/geometry_type/trigonal_bipyramidal.xyz +++ b/tests/testresources/inputs/geometry_type/trigonal_bipyramidal.xyz @@ -1,18 +1,18 @@ 16 -07/15/2022 21:19, XYZ structure generated by mol3D Class, molSimplify +04/30/2024 13:23, XYZ structure generated by mol3D Class, molSimplify Fe 0.000000 0.000000 0.000000 -O -2.120000 -0.000421 -0.000881 -H -2.730515 0.026086 0.779322 -H -2.727124 -0.028395 -0.783794 -O 2.120000 -0.000421 -0.000881 -H 2.728823 -0.572514 -0.533990 -H 2.728652 0.572801 0.531233 -O 0.000000 1.629945 -1.355611 -H 0.153088 1.607265 -2.334413 -H -0.153473 2.588291 -1.155199 -O 0.000000 0.352035 2.090567 -H 0.153533 1.209399 2.563291 -H -0.153710 -0.303124 2.818062 -O 0.000000 -1.984703 -0.745222 -H 0.126817 -2.825828 -0.236778 -H -0.126096 -2.282875 -1.681767 +O 2.120000 0.000000 0.000000 +H 2.726198 0.027327 -0.782547 +H 2.726184 -0.027327 0.782547 +O -1.060000 1.835974 0.000000 +H -1.410402 2.333646 -0.781116 +H -1.315789 2.388255 0.781116 +O -1.060000 -1.835974 0.000000 +H -1.292217 -2.401880 -0.778734 +H -1.433975 -2.320020 0.778734 +O 0.000000 0.000000 2.120000 +H -0.293326 -0.726007 2.726198 +H 0.293326 0.726007 2.726184 +O 0.000000 -0.000000 -2.120000 +H 0.222961 0.750609 -2.726198 +H -0.222961 -0.750609 -2.726184 diff --git a/tests/testresources/inputs/geometry_type/trigonal_planar.xyz b/tests/testresources/inputs/geometry_type/trigonal_planar.xyz index 523e46cd..d57b0608 100644 --- a/tests/testresources/inputs/geometry_type/trigonal_planar.xyz +++ b/tests/testresources/inputs/geometry_type/trigonal_planar.xyz @@ -1,12 +1,12 @@ 10 -07/15/2022 21:12, XYZ structure generated by mol3D Class, molSimplify +04/30/2024 13:14, XYZ structure generated by mol3D Class, molSimplify Fe 0.000000 0.000000 0.000000 -O -1.957585 0.492141 0.648119 -H -2.404883 0.225459 1.491371 -H -2.634735 1.041484 0.176951 -O 0.951527 -1.746264 0.734547 -H 1.923472 -1.917547 0.825504 -H 0.526130 -2.578001 1.065509 -O 0.965381 1.223677 -1.437030 -H 0.569473 1.950860 -1.981760 -H 1.915622 1.199311 -1.717647 +O 2.120000 0.000000 0.000000 +H 2.726198 0.027327 -0.782547 +H 2.726184 -0.027327 0.782547 +O -1.060000 1.835974 0.000000 +H -1.410402 2.333646 -0.781116 +H -1.315789 2.388255 0.781116 +O -1.060000 -1.835974 0.000000 +H -1.292217 -2.401880 -0.778734 +H -1.433975 -2.320020 0.778734