From 6601dc5b7328ce2330706832e60ee97822fdd630 Mon Sep 17 00:00:00 2001 From: aarongarrison Date: Tue, 26 Mar 2024 09:11:45 -0400 Subject: [PATCH 01/12] Adding geometry classification based on shape matching. --- molSimplify/Classes/globalvars.py | 31 ++++++ molSimplify/Classes/mol3D.py | 158 +++++++++++++++++++++++++++++- 2 files changed, 188 insertions(+), 1 deletion(-) diff --git a/molSimplify/Classes/globalvars.py b/molSimplify/Classes/globalvars.py index 54e45098..44efc670 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,26 @@ [[70.5, 70.5, 67.7, 67.7, 120, 120, 135.5, 135.5]] } +#need these 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)]) +} + # Module for running bash commands # @param cmd String containing command to be run # @return bash output string @@ -720,6 +741,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..c31d6f5d 100644 --- a/molSimplify/Classes/mol3D.py +++ b/molSimplify/Classes/mol3D.py @@ -25,7 +25,7 @@ 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 try: import PyQt5 # noqa: F401 @@ -5683,6 +5683,162 @@ def get_geometry_type(self, dict_check=False, angle_ref=False, } return results + def get_geometry_type_distance(self, max_dev=1e6, + flag_catoms=False, catoms_arr=None, debug=False, + skip=False, transition_metals_only=False): + """ + 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. + 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 + Contains the classified geometry and the RMSD from an ideal structure. + Summary contains the RMSD 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: + 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 list(all_geometries.keys()): + #should we indicate somehow that these are unknown due to a different coordination number? + results = { + "geometry": "unknown", + "rmsd": False, + "summary": {}, + "hapticity": hapt, + } + return results + + possible_geometries = all_geometries[num_coord] + + for geotype in possible_geometries: + rmsd_calc = self.dev_from_ideal_geometry(all_polyhedra[geotype]) + if debug: + print("Geocheck assigned catoms: ", catoms_assigned, + [first_shell.getAtom(ind).symbol() for ind in catoms_assigned]) + summary.update({geotype: rmsd_calc}) + + current_rmsd, geometry = max_dev, None + for geotype in summary: + if summary[geotype] < current_rmsd: + current_rmsd = summary[geotype] + geometry = geotype + results = { + "geometry": geometry, + "rmsd": current_rmsd, + "summary": summary, + "hapticity": hapt, + } + return results + + def dev_from_ideal_geometry(self, ideal_polyhedron, max_dev=1e6): + """ + 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 + max_dev : float, optional + Maximum RMSD allowed between a structure and an ideal geometry before it is classified as unknown. Default is 1e6. + + Returns + ------- + rmsd: float + Minimum root mean square distance between the fed geometry and the ideal polyhedron + """ + + #get the average bond distance for the first coordination shell of the provided geometry + 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.') + fcs_indices = self.get_fcs() + if len(fcs_indices) - 1 != len(ideal_polyhedron): + raise ValueError('The coordination number differs between the two provided structures.') + metal_atom = self.getAtoms()[metal_idx[0]] + fcs_atoms = [self.getAtoms()[i] for i in fcs_indices] + distances = [] + positions = np.zeros([len(fcs_indices), 3]) + for idx, atom in enumerate(fcs_atoms): + if atom is not metal_atom: + distance = atom.distance(metal_atom) + distances.append(distance) + positions[idx, :] = atom.coords() + mean_dist = np.array(distances).mean() + + #scale ideal geometry to have same average distance + scaled_polyhedron = ideal_polyhedron * mean_dist + + def permutations(list): + 'Returns all possible permutations of a list.' + if len(list) == 0: + return [] + elif len(list) == 1: + return [list] + l = [] + for i in range(len(list)): + m = list[i] + remaining = list[:i] + list[i+1:] + for p in permutations(remaining): + l.append([m] + p) + return l + + current_min = max_dev + orders = permutations(list(np.arange(1, len(ideal_polyhedron)+1))) + + ideal_positions = np.zeros([len(fcs_indices), 3]) + for order in orders: + for i in range(len(order)): + ideal_positions[i+1, :] = scaled_polyhedron[order[i]-1] + rmsd_calc = kabsch_rmsd(ideal_positions, positions) + if rmsd_calc < current_min: + current_min = rmsd_calc + + return current_min + + + #TODO: build a mol3D object from a list of 3D coordinates with the metal at (0, 0, 0) + #TODO: calculate the average bond distance for the provided mol3D + #TODO: scale the ideal polyhedra by the average bond distance + #TODO: calculate the RMSD, return the minimum + 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): From bfef0d5ef4b38f5863ee2c0fffce1bb8fbf8a52b Mon Sep 17 00:00:00 2001 From: aarongarrison Date: Tue, 2 Apr 2024 11:00:07 -0400 Subject: [PATCH 02/12] Adding translation, cleaning up comments --- molSimplify/Classes/globalvars.py | 2 +- molSimplify/Classes/mol3D.py | 22 ++++++++-------------- 2 files changed, 9 insertions(+), 15 deletions(-) diff --git a/molSimplify/Classes/globalvars.py b/molSimplify/Classes/globalvars.py index 44efc670..0e4678c0 100644 --- a/molSimplify/Classes/globalvars.py +++ b/molSimplify/Classes/globalvars.py @@ -464,7 +464,7 @@ [[70.5, 70.5, 67.7, 67.7, 120, 120, 135.5, 135.5]] } -#need these to be distance 1 from the origin +#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)]), diff --git a/molSimplify/Classes/mol3D.py b/molSimplify/Classes/mol3D.py index c31d6f5d..5192e832 100644 --- a/molSimplify/Classes/mol3D.py +++ b/molSimplify/Classes/mol3D.py @@ -5748,8 +5748,8 @@ def get_geometry_type_distance(self, max_dev=1e6, for geotype in possible_geometries: rmsd_calc = self.dev_from_ideal_geometry(all_polyhedra[geotype]) if debug: - print("Geocheck assigned catoms: ", catoms_assigned, - [first_shell.getAtom(ind).symbol() for ind in catoms_assigned]) + pass + #not sure what to include here summary.update({geotype: rmsd_calc}) current_rmsd, geometry = max_dev, None @@ -5765,7 +5765,7 @@ def get_geometry_type_distance(self, max_dev=1e6, } return results - def dev_from_ideal_geometry(self, ideal_polyhedron, max_dev=1e6): + def dev_from_ideal_geometry(self, ideal_polyhedron): """ 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. @@ -5774,8 +5774,6 @@ def dev_from_ideal_geometry(self, ideal_polyhedron, max_dev=1e6): ---------- ideal_polyhedron: np.array of 3-tuples of coordinates Reference list of points for an ideal geometry - max_dev : float, optional - Maximum RMSD allowed between a structure and an ideal geometry before it is classified as unknown. Default is 1e6. Returns ------- @@ -5783,7 +5781,6 @@ def dev_from_ideal_geometry(self, ideal_polyhedron, max_dev=1e6): Minimum root mean square distance between the fed geometry and the ideal polyhedron """ - #get the average bond distance for the first coordination shell of the provided geometry metal_idx = self.findMetal() if len(metal_idx) == 0: raise ValueError('No metal centers exist in this complex.') @@ -5794,13 +5791,14 @@ def dev_from_ideal_geometry(self, ideal_polyhedron, max_dev=1e6): raise ValueError('The coordination number differs between the two provided structures.') metal_atom = self.getAtoms()[metal_idx[0]] fcs_atoms = [self.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): if atom is not metal_atom: distance = atom.distance(metal_atom) distances.append(distance) - positions[idx, :] = atom.coords() + positions[idx, :] = atom.coords() - metal_atom.coords() #shift so the metal is at (0, 0, 0) mean_dist = np.array(distances).mean() #scale ideal geometry to have same average distance @@ -5820,9 +5818,10 @@ def permutations(list): l.append([m] + p) return l - current_min = max_dev + current_min = np.inf orders = permutations(list(np.arange(1, len(ideal_polyhedron)+1))) + #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)): @@ -5831,13 +5830,8 @@ def permutations(list): if rmsd_calc < current_min: current_min = rmsd_calc + #return minimum RMSD return current_min - - - #TODO: build a mol3D object from a list of 3D coordinates with the metal at (0, 0, 0) - #TODO: calculate the average bond distance for the provided mol3D - #TODO: scale the ideal polyhedra by the average bond distance - #TODO: calculate the RMSD, return the minimum def get_features(self, lac=True, force_generate=False, eq_sym=False, use_dist=False, NumB=False, Gval=False, size_normalize=False, From b576b861efbf33ae2fcf561958604833c2846877 Mon Sep 17 00:00:00 2001 From: jwtoney Date: Wed, 3 Apr 2024 16:08:16 -0400 Subject: [PATCH 03/12] fixed geometry classification --- molSimplify/Classes/globalvars.py | 4 +++- molSimplify/Classes/mol3D.py | 22 +++++++++------------- 2 files changed, 12 insertions(+), 14 deletions(-) diff --git a/molSimplify/Classes/globalvars.py b/molSimplify/Classes/globalvars.py index 0e4678c0..6ef27f9c 100644 --- a/molSimplify/Classes/globalvars.py +++ b/molSimplify/Classes/globalvars.py @@ -481,7 +481,9 @@ "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)]) + "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 diff --git a/molSimplify/Classes/mol3D.py b/molSimplify/Classes/mol3D.py index 5192e832..27f69510 100644 --- a/molSimplify/Classes/mol3D.py +++ b/molSimplify/Classes/mol3D.py @@ -1979,7 +1979,7 @@ def getBondedAtomsOct(self, ind, CN=6, debug=False, flag_loose=False, atom_speci 'Error, mol3D could not understand connectivity in mol') return nats - def getBondedAtomsSmart(self, idx, oct=True, strict_cutoff=False, catom_list=None): + def getBondedAtomsSmart(self, idx, oct=False, strict_cutoff=False, catom_list=None): """ Get the atoms bonded with the atom specified with the given index, using the molecular graph. Creates graph if it does not exist. @@ -4938,8 +4938,8 @@ 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: - _, catoms = self.oct_comp(debug=False) + # if len(catoms) > 6: + # _, catoms = self.oct_comp(debug=False) fcs = [metalind] + catoms return fcs @@ -5684,7 +5684,7 @@ def get_geometry_type(self, dict_check=False, angle_ref=False, return results def get_geometry_type_distance(self, max_dev=1e6, - flag_catoms=False, catoms_arr=None, debug=False, + flag_catoms=False, catoms_arr=None, skip=False, transition_metals_only=False): """ Get the type of the geometry (available options in globalvars all_geometries). @@ -5701,8 +5701,6 @@ def get_geometry_type_distance(self, max_dev=1e6, 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 @@ -5729,7 +5727,6 @@ def get_geometry_type_distance(self, max_dev=1e6, 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.") @@ -5747,12 +5744,9 @@ def get_geometry_type_distance(self, max_dev=1e6, for geotype in possible_geometries: rmsd_calc = self.dev_from_ideal_geometry(all_polyhedra[geotype]) - if debug: - pass - #not sure what to include here summary.update({geotype: rmsd_calc}) - current_rmsd, geometry = max_dev, None + current_rmsd, geometry = max_dev, "unknown" for geotype in summary: if summary[geotype] < current_rmsd: current_rmsd = summary[geotype] @@ -5786,7 +5780,9 @@ def dev_from_ideal_geometry(self, ideal_polyhedron): raise ValueError('No metal centers exist in this complex.') elif len(metal_idx) != 1: raise ValueError('Multimetal complexes are not yet handled.') - fcs_indices = self.get_fcs() + temp_mol = self.get_first_shell() + fcs_indices = temp_mol[0].get_fcs() + if len(fcs_indices) - 1 != len(ideal_polyhedron): raise ValueError('The coordination number differs between the two provided structures.') metal_atom = self.getAtoms()[metal_idx[0]] @@ -5798,7 +5794,7 @@ def dev_from_ideal_geometry(self, ideal_polyhedron): if atom is not metal_atom: distance = atom.distance(metal_atom) distances.append(distance) - positions[idx, :] = atom.coords() - metal_atom.coords() #shift so the metal is at (0, 0, 0) + positions[idx, :] = np.array(atom.coords()) - np.array(metal_atom.coords()) #shift so the metal is at (0, 0, 0) mean_dist = np.array(distances).mean() #scale ideal geometry to have same average distance From d36f4c75d74b41b0d0a6105fed8d633f2b2a799e Mon Sep 17 00:00:00 2001 From: jwtoney Date: Thu, 4 Apr 2024 17:21:54 -0400 Subject: [PATCH 04/12] debugging get_geometry_distance() --- molSimplify/Classes/mol3D.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/molSimplify/Classes/mol3D.py b/molSimplify/Classes/mol3D.py index 27f69510..8810e3b3 100644 --- a/molSimplify/Classes/mol3D.py +++ b/molSimplify/Classes/mol3D.py @@ -5780,21 +5780,27 @@ def dev_from_ideal_geometry(self, ideal_polyhedron): 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() - fcs_indices = temp_mol[0].get_fcs() + temp_mol = self.get_first_shell()[0] + fcs_indices = temp_mol.get_fcs() + # remove metal index from first coordination shell + fcs_indices.pop(temp_mol.findMetal()[0]) - if len(fcs_indices) - 1 != len(ideal_polyhedron): + if len(fcs_indices) != len(ideal_polyhedron): raise ValueError('The coordination number differs between the two provided structures.') metal_atom = self.getAtoms()[metal_idx[0]] + print(metal_atom.coords()) fcs_atoms = [self.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 atom in fcs_atoms: + print(atom.coords()) for idx, atom in enumerate(fcs_atoms): if atom is not metal_atom: distance = atom.distance(metal_atom) distances.append(distance) positions[idx, :] = np.array(atom.coords()) - np.array(metal_atom.coords()) #shift so the metal is at (0, 0, 0) + print(positions) mean_dist = np.array(distances).mean() #scale ideal geometry to have same average distance @@ -5815,14 +5821,14 @@ def permutations(list): return l current_min = np.inf - orders = permutations(list(np.arange(1, len(ideal_polyhedron)+1))) + orders = permutations(list(np.arange(1, len(ideal_polyhedron)))) #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)): - ideal_positions[i+1, :] = scaled_polyhedron[order[i]-1] - rmsd_calc = kabsch_rmsd(ideal_positions, positions) + ideal_positions[i, :] = scaled_polyhedron[order[i]-1] + rmsd_calc = kabsch_rmsd(ideal_positions, positions, translate=True) if rmsd_calc < current_min: current_min = rmsd_calc From 63951334195ad79ec104d1ab012c58bc206786cb Mon Sep 17 00:00:00 2001 From: aarongarrison Date: Thu, 4 Apr 2024 19:49:59 -0400 Subject: [PATCH 05/12] Fixing an indexing error introduced in the edge/sandwich fix --- molSimplify/Classes/mol3D.py | 35 +++++++++++++++++------------------ 1 file changed, 17 insertions(+), 18 deletions(-) diff --git a/molSimplify/Classes/mol3D.py b/molSimplify/Classes/mol3D.py index 8810e3b3..ffdbd4c2 100644 --- a/molSimplify/Classes/mol3D.py +++ b/molSimplify/Classes/mol3D.py @@ -5783,28 +5783,25 @@ def dev_from_ideal_geometry(self, ideal_polyhedron): temp_mol = self.get_first_shell()[0] fcs_indices = temp_mol.get_fcs() # remove metal index from first coordination shell - fcs_indices.pop(temp_mol.findMetal()[0]) + 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.') - metal_atom = self.getAtoms()[metal_idx[0]] - print(metal_atom.coords()) - fcs_atoms = [self.getAtoms()[i] for i in fcs_indices] + + #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 atom in fcs_atoms: - print(atom.coords()) for idx, atom in enumerate(fcs_atoms): - if atom is not metal_atom: - distance = atom.distance(metal_atom) - distances.append(distance) - positions[idx, :] = np.array(atom.coords()) - np.array(metal_atom.coords()) #shift so the metal is at (0, 0, 0) - print(positions) - mean_dist = np.array(distances).mean() - - #scale ideal geometry to have same average distance - scaled_polyhedron = ideal_polyhedron * mean_dist + distance = atom.distance(metal_atom) + distances.append(distance) + positions[idx, :] = np.array(atom.coords()) - np.array(metal_atom.coords()) #shift so the metal is at (0, 0, 0) + + #make it so the ideal polyhedron has same average bond distance as mol + scaled_polyhedron = ideal_polyhedron * np.mean(np.array(distances)) def permutations(list): 'Returns all possible permutations of a list.' @@ -5821,14 +5818,16 @@ def permutations(list): return l current_min = np.inf - orders = permutations(list(np.arange(1, len(ideal_polyhedron)))) + orders = permutations(list(range(len(ideal_polyhedron)))) #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)): - ideal_positions[i, :] = scaled_polyhedron[order[i]-1] - rmsd_calc = kabsch_rmsd(ideal_positions, positions, translate=True) + ideal_positions[i, :] = scaled_polyhedron[order[i]] + #if you wanted to let each ligand scale its bond 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 From 391d51c7e659de482f1b5f50546c67250530488a Mon Sep 17 00:00:00 2001 From: aarongarrison Date: Fri, 5 Apr 2024 14:02:48 -0400 Subject: [PATCH 06/12] Add a tiebreaker metric and allowed ligand lengths to scale independently --- molSimplify/Classes/mol3D.py | 60 +++++++++++++++++++++++++----------- 1 file changed, 42 insertions(+), 18 deletions(-) diff --git a/molSimplify/Classes/mol3D.py b/molSimplify/Classes/mol3D.py index ffdbd4c2..25479d45 100644 --- a/molSimplify/Classes/mol3D.py +++ b/molSimplify/Classes/mol3D.py @@ -25,7 +25,7 @@ from molSimplify.Scripts.geometry import (distance, connectivity_match, vecangle, rotation_params, rotate_around_axis) -from molSimplify.Scripts.rmsd import rigorous_rmsd, kabsch_rmsd +from molSimplify.Scripts.rmsd import rigorous_rmsd, kabsch_rmsd, kabsch_rotate try: import PyQt5 # noqa: F401 @@ -5683,7 +5683,7 @@ def get_geometry_type(self, dict_check=False, angle_ref=False, } return results - def get_geometry_type_distance(self, max_dev=1e6, + 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): """ @@ -5696,6 +5696,8 @@ def get_geometry_type_distance(self, max_dev=1e6, ---------- 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 @@ -5710,7 +5712,7 @@ def get_geometry_type_distance(self, max_dev=1e6, ------- results : dictionary Contains the classified geometry and the RMSD from an ideal structure. - Summary contains the RMSD for all considered geometry types. + Summary contains a list of the RMSD and the maximum single-atom deviation for all considered geometry types. """ @@ -5741,21 +5743,34 @@ def get_geometry_type_distance(self, max_dev=1e6, 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 = self.dev_from_ideal_geometry(all_polyhedra[geotype]) - summary.update({geotype: rmsd_calc}) + 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 summary[geotype] < current_rmsd: - current_rmsd = summary[geotype] - geometry = geotype + #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 @@ -5773,6 +5788,8 @@ def dev_from_ideal_geometry(self, ideal_polyhedron): ------- 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() @@ -5799,9 +5816,6 @@ def dev_from_ideal_geometry(self, ideal_polyhedron): distance = atom.distance(metal_atom) distances.append(distance) positions[idx, :] = np.array(atom.coords()) - np.array(metal_atom.coords()) #shift so the metal is at (0, 0, 0) - - #make it so the ideal polyhedron has same average bond distance as mol - scaled_polyhedron = ideal_polyhedron * np.mean(np.array(distances)) def permutations(list): 'Returns all possible permutations of a list.' @@ -5819,20 +5833,30 @@ def permutations(list): current_min = np.inf orders = permutations(list(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)): - ideal_positions[i, :] = scaled_polyhedron[order[i]] - #if you wanted to let each ligand scale its bond length independently, uncomment the following - #ideal_positions[i, :] = ideal_polyhedron[order[i]] * distances[i] + #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 - - #return minimum RMSD - return current_min + #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, From 51dd7b7e71dd5a5fb79e9018a82435a6bc2abf4a Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 9 Apr 2024 20:23:45 +0000 Subject: [PATCH 07/12] [pre-commit.ci] auto fixes from pre-commit hooks --- molSimplify/Classes/mol3D.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/molSimplify/Classes/mol3D.py b/molSimplify/Classes/mol3D.py index 25479d45..dc724445 100644 --- a/molSimplify/Classes/mol3D.py +++ b/molSimplify/Classes/mol3D.py @@ -5764,7 +5764,7 @@ def get_geometry_type_distance(self, max_dev=1e6, close_dev=1e-2, else: current_rmsd = summary[geotype][0] geometry = geotype - + results = { "geometry": geometry, "rmsd": current_rmsd, @@ -5816,7 +5816,7 @@ def dev_from_ideal_geometry(self, ideal_polyhedron): distance = atom.distance(metal_atom) distances.append(distance) positions[idx, :] = np.array(atom.coords()) - np.array(metal_atom.coords()) #shift so the metal is at (0, 0, 0) - + def permutations(list): 'Returns all possible permutations of a list.' if len(list) == 0: @@ -5837,7 +5837,7 @@ def permutations(list): #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: From 79272c199d4c67ccd4bcb43c2e11b57195d113c8 Mon Sep 17 00:00:00 2001 From: aarongarrison Date: Tue, 16 Apr 2024 16:25:09 -0400 Subject: [PATCH 08/12] Returning GetBondedAtomsSmart default to octahedral --- molSimplify/Classes/mol3D.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/molSimplify/Classes/mol3D.py b/molSimplify/Classes/mol3D.py index 25479d45..a1b8ed2a 100644 --- a/molSimplify/Classes/mol3D.py +++ b/molSimplify/Classes/mol3D.py @@ -1979,7 +1979,7 @@ def getBondedAtomsOct(self, ind, CN=6, debug=False, flag_loose=False, atom_speci 'Error, mol3D could not understand connectivity in mol') return nats - def getBondedAtomsSmart(self, idx, oct=False, strict_cutoff=False, catom_list=None): + def getBondedAtomsSmart(self, idx, oct=True, strict_cutoff=False, catom_list=None): """ Get the atoms bonded with the atom specified with the given index, using the molecular graph. Creates graph if it does not exist. @@ -5726,7 +5726,8 @@ def get_geometry_type_distance(self, max_dev=1e6, close_dev=1e-2, 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])) + #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: From 46abd16ff414df11fff9c5b18b6e762e01a38483 Mon Sep 17 00:00:00 2001 From: aarongarrison Date: Tue, 16 Apr 2024 16:59:19 -0400 Subject: [PATCH 09/12] Reverted behavior of get_fcs back to default. Added in a keyword argument to allow for higher than 6-coordinated structures. --- molSimplify/Classes/mol3D.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/molSimplify/Classes/mol3D.py b/molSimplify/Classes/mol3D.py index 2b62fc2f..76640465 100644 --- a/molSimplify/Classes/mol3D.py +++ b/molSimplify/Classes/mol3D.py @@ -4918,7 +4918,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 +4928,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,8 +4940,8 @@ 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: - # _, catoms = self.oct_comp(debug=False) + if max6 and len(catoms) > 6: + _, catoms = self.oct_comp(debug=False) fcs = [metalind] + catoms return fcs @@ -5727,7 +5729,7 @@ def get_geometry_type_distance(self, max_dev=1e6, close_dev=1e-2, 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) + 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: @@ -5799,7 +5801,7 @@ def dev_from_ideal_geometry(self, ideal_polyhedron): 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() + fcs_indices = temp_mol.get_fcs(max6=False) # remove metal index from first coordination shell fcs_indices.remove(temp_mol.findMetal()[0]) From 60122e13281a0cc9a8da6f8ff0f632b35248c4c6 Mon Sep 17 00:00:00 2001 From: aarongarrison Date: Fri, 26 Apr 2024 09:43:53 -0400 Subject: [PATCH 10/12] Updating format of unknown geometry summaries to match other summaries --- molSimplify/Classes/mol3D.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/molSimplify/Classes/mol3D.py b/molSimplify/Classes/mol3D.py index 76640465..2b49828d 100644 --- a/molSimplify/Classes/mol3D.py +++ b/molSimplify/Classes/mol3D.py @@ -5739,9 +5739,10 @@ def get_geometry_type_distance(self, max_dev=1e6, close_dev=1e-2, #should we indicate somehow that these are unknown due to a different coordination number? results = { "geometry": "unknown", - "rmsd": False, + "rmsd": np.NAN, "summary": {}, "hapticity": hapt, + "close_rmsds": False } return results From 9378d277d32025c32bad1eb7d0e497bbb02bad7b Mon Sep 17 00:00:00 2001 From: Ralf Meyer Date: Tue, 30 Apr 2024 13:40:43 -0400 Subject: [PATCH 11/12] Add typing and switch to itertools.permutation --- molSimplify/Classes/mol3D.py | 67 +++++++++++++++--------------------- 1 file changed, 28 insertions(+), 39 deletions(-) diff --git a/molSimplify/Classes/mol3D.py b/molSimplify/Classes/mol3D.py index 2b49828d..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 @@ -26,6 +26,7 @@ vecangle, rotation_params, rotate_around_axis) from molSimplify.Scripts.rmsd import rigorous_rmsd, kabsch_rmsd, kabsch_rotate +from itertools import permutations try: import PyQt5 # noqa: F401 @@ -5685,9 +5686,10 @@ 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): + 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). @@ -5728,7 +5730,7 @@ def get_geometry_type_distance(self, max_dev=1e6, close_dev=1e-2, 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 + # 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.') @@ -5736,19 +5738,19 @@ def get_geometry_type_distance(self, max_dev=1e6, close_dev=1e-2, 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? + # 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 + "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 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]}) @@ -5756,13 +5758,13 @@ def get_geometry_type_distance(self, max_dev=1e6, close_dev=1e-2, 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 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 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 + # classify based on largest singular deviation current_rmsd = summary[geotype][0] geometry = geotype else: @@ -5778,7 +5780,7 @@ def get_geometry_type_distance(self, max_dev=1e6, close_dev=1e-2, } return results - def dev_from_ideal_geometry(self, ideal_polyhedron): + 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. @@ -5809,57 +5811,44 @@ def dev_from_ideal_geometry(self, ideal_polyhedron): 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 + # 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 + # 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) - positions[idx, :] = np.array(atom.coords()) - np.array(metal_atom.coords()) #shift so the metal is at (0, 0, 0) - - def permutations(list): - 'Returns all possible permutations of a list.' - if len(list) == 0: - return [] - elif len(list) == 1: - return [list] - l = [] - for i in range(len(list)): - m = list[i] - remaining = list[:i] + list[i+1:] - for p in permutations(remaining): - l.append([m] + p) - return l + # 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(list(range(len(ideal_polyhedron)))) + 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)) + # 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 + # 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 + # 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 + # 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 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, From ae35599d5b65706aac4da2e9e176740a8d0f4e68 Mon Sep 17 00:00:00 2001 From: Ralf Meyer Date: Tue, 30 Apr 2024 13:41:12 -0400 Subject: [PATCH 12/12] Add test cases and fix reference geometries --- molSimplify/Data/tbp.dat | 12 ++-- molSimplify/Data/tpl.dat | 6 +- tests/test_mol3D.py | 66 +++++++++++++++++++ .../inputs/geometry_type/seesaw.xyz | 26 ++++---- .../geometry_type/trigonal_bipyramidal.xyz | 32 ++++----- .../inputs/geometry_type/trigonal_planar.xyz | 20 +++--- 6 files changed, 114 insertions(+), 48 deletions(-) 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