From f4ba1115bb9ceb874dc8a9645227a7acab941c62 Mon Sep 17 00:00:00 2001 From: jwtoney Date: Tue, 10 Dec 2024 19:45:43 -0500 Subject: [PATCH 1/3] lots of code to identity and flip TMCs to new symmetry groups --- molSimplify/Classes/mol3D.py | 488 +++++++++++++++++++++++++++++++++-- 1 file changed, 463 insertions(+), 25 deletions(-) diff --git a/molSimplify/Classes/mol3D.py b/molSimplify/Classes/mol3D.py index 4f7259fb..4d69f1da 100644 --- a/molSimplify/Classes/mol3D.py +++ b/molSimplify/Classes/mol3D.py @@ -6084,24 +6084,32 @@ def dev_from_ideal_geometry(self, ideal_polyhedron: np.ndarray) -> Tuple[float, # return minimum RMSD, maximum pairwise distance in that structure return current_min, max_dist - def get_symmetry(tmc_mol, verbose=True, max_allowed_dev=30): + def get_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, details=False): """ - Classify transition metal complexes (TMCs) according to symmetry group + Classify octahedral transition metal complexes (TMCs) according to symmetry Parameters ---------- - tmc_mol : mol3D - mol3D instance of TMC with unknown symmetry group. + tmc_mol: mol3D + mol3D instance of TMC with unknown symmetry verbose: bool - Flag for returning warning when TMC exhibits high deviation from closest symmetry group. Default is True. + Flag for returning warning when TMC exhibits high deviation from closest symmetry + Default=True max_allowed_dev: float - Maximum allowed deviation before warning is triggered (degrees). Default is 30. + Maximum allowed deviation before warning is triggered (degrees) + Default=30 + details: bool + Flag for returning detailed lists of unique ligands and coordinating atoms, intended for use with flip_symmetry() + Default=False Returns ------- - symmetry_dict : dictionary - Dictionary storing assigned symmetry class and devations from all possible symmetry classes. + symmetry_dict: dict + Dictionary storing assigned symmetry class and deviations from all possible symmetry classes + detailed_dict: dict + Dictionary storing indices of metal, ligand atoms, and ligand coordinating atoms, only returned if details=True """ + from molSimplify.Classes.ligand import ligand_breakdown from molSimplify.Classes.mol2D import Mol2D metal_idx = tmc_mol.findMetal() @@ -6126,7 +6134,8 @@ def get_symmetry(tmc_mol, verbose=True, max_allowed_dev=30): ligand_mol.deleteatoms(Alist=[i for i in tmc_atoms if i not in ligand]) ligand_hashes.append(Mol2D().from_mol3d(mol3d=ligand_mol).graph_hash()) - # determine number of unique ligands + # determine number of unique ligands, sorted in ascending order + # sort is stable, so in the case of ties, original order is preserved unique_ligands = dict(sorted(Counter(ligand_hashes).items(), key=lambda x: x[1])) if len(unique_ligands) > 3: raise ValueError('Function only supported for TMCs with up to 3 unique ligands.') @@ -6137,18 +6146,25 @@ def get_symmetry(tmc_mol, verbose=True, max_allowed_dev=30): symmetry_dict = {'symmetry': symmetry} # two unique ligands: consider cis, trans, fac, mer, and monoheteroleptic symmetry groups - if len(unique_ligands) == 2: + elif len(unique_ligands) == 2: # calculate ratio between ligands - unique_ligand_ratio = list(unique_ligands.values())[1] / list(unique_ligands.values())[0] - lig2_indices = [index for index, value in enumerate(ligand_hashes) if - value == list(unique_ligands.keys())[0]] + unique_ligand_ratios = list(unique_ligands.values())[1] / list(unique_ligands.values())[0] + lig2_indices = [index for index, value in enumerate(ligand_hashes) if value == list(unique_ligands.keys())[0]] + lig2_atoms = [lig_list[idx] for idx in lig2_indices] lig2_catoms = [lig_catoms[idx][0] for idx in lig2_indices] - if unique_ligand_ratio == 5: + lig1_indices = [index for index, value in enumerate(ligand_hashes) if value == list(unique_ligands.keys())[1]] + lig1_atoms = [lig_list[idx] for idx in lig1_indices] + lig1_catoms = [lig_catoms[idx][0] for idx in lig1_indices] + + lig3_atoms = False + lig3_catoms = False + + if unique_ligand_ratios == 5: symmetry = 'monoheteroleptic' symmetry_dict = {'symmetry': symmetry} - elif unique_ligand_ratio == 2: + elif unique_ligand_ratios == 2: # find angle between coordinating atoms of less represented (i.e., minority) ligand and metal lig2_angle = tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig2_catoms[1]) # classify complex as cis or trans based on deviation from ideal angle @@ -6162,7 +6178,7 @@ def get_symmetry(tmc_mol, verbose=True, max_allowed_dev=30): print('Warning: high deviation from ideal symmetry, manual inspection recommended') symmetry_dict = {'symmetry': symmetry, 'cis_dev': cis_dev, 'trans_dev': trans_dev} - elif unique_ligand_ratio == 1: + elif unique_ligand_ratios == 1: # find angle between coordinating atoms of first ligand and metal lig2_angles = np.array((tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig2_catoms[1]), tmc_mol.getAngle(idx0=lig2_catoms[1], idx1=metal_idx, idx2=lig2_catoms[2]), @@ -6183,8 +6199,8 @@ def get_symmetry(tmc_mol, verbose=True, max_allowed_dev=30): # three unique ligands: consider cis asymmetric (CA), double cis symmetric (DCS), trans asymmetric (TA), # double trans symmetric (DTS), equatorial asymmetric (EA), fac asymmetric (FA), mer asymmetric trans (MAT), # and mer asymmetric cis (MAC) symmetry groups - if len(unique_ligands) == 3: - # calculate ratio between ligands + elif len(unique_ligands) == 3: + # calculate ratio between ligands sorted in ascending order (L1 = unique_ligands[2] = most abundant) # ratios stored as follows: L1:L2, L2:L3, L1:L3 unique_ligand_ratios = [list(unique_ligands.values())[2] / list(unique_ligands.values())[1], list(unique_ligands.values())[1] / list(unique_ligands.values())[0], @@ -6192,14 +6208,17 @@ def get_symmetry(tmc_mol, verbose=True, max_allowed_dev=30): lig1_indices = [index for index, value in enumerate(ligand_hashes) if value == list(unique_ligands.keys())[2]] + lig1_atoms = [lig_list[idx] for idx in lig1_indices] lig1_catoms = [lig_catoms[idx][0] for idx in lig1_indices] lig2_indices = [index for index, value in enumerate(ligand_hashes) if value == list(unique_ligands.keys())[1]] + lig2_atoms = [lig_list[idx] for idx in lig2_indices] lig2_catoms = [lig_catoms[idx][0] for idx in lig2_indices] lig3_indices = [index for index, value in enumerate(ligand_hashes) if value == list(unique_ligands.keys())[0]] + lig3_atoms = [lig_list[idx] for idx in lig3_indices] lig3_catoms = [lig_catoms[idx][0] for idx in lig3_indices] if unique_ligand_ratios == [4, 1, 4]: @@ -6224,19 +6243,19 @@ def get_symmetry(tmc_mol, verbose=True, max_allowed_dev=30): lig_angles.sort() # classify complex as DCS, DCT, or EA based on deviation from ideal angle DCS_dev_avg = np.average(np.abs(lig_angles - 90)) - DCT_dev_avg = np.average(np.abs(lig_angles - 180)) + DTS_dev_avg = np.average(np.abs(lig_angles - 180)) EA_dev_avg = np.average( np.abs(np.concatenate((np.array(lig_angles)[0:2] - 90, np.array([lig_angles[2] - 180]))))) - if min(DCS_dev_avg, DCT_dev_avg, EA_dev_avg) == DCS_dev_avg: + if min(DCS_dev_avg, DTS_dev_avg, EA_dev_avg) == DCS_dev_avg: symmetry = 'double cis symmetric' - elif min(DCS_dev_avg, DCT_dev_avg, EA_dev_avg) == DCT_dev_avg: + elif min(DCS_dev_avg, DTS_dev_avg, EA_dev_avg) == DTS_dev_avg: symmetry = 'double trans symmetric' - elif min(DCS_dev_avg, DCT_dev_avg, EA_dev_avg) == EA_dev_avg: + elif min(DCS_dev_avg, DTS_dev_avg, EA_dev_avg) == EA_dev_avg: symmetry = 'equatorial asymmetric' - if verbose and min(DCS_dev_avg, DCT_dev_avg, EA_dev_avg) > max_allowed_dev: + if verbose and min(DCS_dev_avg, DTS_dev_avg, EA_dev_avg) > max_allowed_dev: print('Warning: high deviation from ideal symmetry, manual inspection recommended') symmetry_dict = {'symmetry': symmetry, 'double_cis_symmetric_dev': DCS_dev_avg, - 'double_trans_symmetric_dev': DCS_dev_avg, 'equatorial_asymmetric_dev': EA_dev_avg} + 'double_trans_symmetric_dev': DTS_dev_avg, 'equatorial_asymmetric_dev': EA_dev_avg} elif unique_ligand_ratios == [3 / 2, 2, 3]: # find all angles between coordinating atoms of all L1 and L2 type ligands and metal @@ -6263,7 +6282,426 @@ def get_symmetry(tmc_mol, verbose=True, max_allowed_dev=30): symmetry_dict = {'symmetry': symmetry, 'fac_asymmetric_dev': FA_dev_avg, 'mer_asymmetric_trans_dev': MAT_dev_avg, 'mer_asymmetric_cis_dev': MAC_dev_avg} - return symmetry_dict + if details: + detailed_dict = {'num_unique_ligands': len(unique_ligands), 'unique_ligand_ratios': unique_ligand_ratios, + 'metal_idx': metal_idx, 'lig1_atoms': lig1_atoms, 'lig1_catoms': lig1_catoms, + 'lig2_atoms': lig2_atoms, 'lig2_catoms': lig2_catoms, + 'lig3_atoms': lig3_atoms, 'lig3_catoms': lig3_catoms} + return symmetry_dict, detailed_dict + else: + return symmetry_dict + + def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=False): + """ + Flip octahedral transition metal complexes (TMCs) to opposite symmetry group. + Example: cis to trans, fac to mer, etc. + + Parameters + ---------- + tmc_mol: mol3D + mol3D instance of TMC + verbose: bool + Flag for returning warning when TMC exhibits high deviation from closest symmetry + Default=True + max_allowed_dev: float + Maximum allowed deviation before warning is triggered (degrees) + Default=30 + target_symmetry: str + Target symmetry for complex to be transformed to, only defined for num_unique_ligands == 3 + Default=False + + Returns + ------- + tmc_mol: mol3D + returns self, a mol3D object with flipped symmetry + """ + + symmetry_dict, detailed_dict, = tmc_mol.get_symmetry(verbose=verbose, + max_allowed_dev=max_allowed_dev, + details=True) + symmetry = symmetry_dict['symmetry'] + + num_unique_ligands = detailed_dict['num_unique_ligands'] + unique_ligand_ratios = detailed_dict['unique_ligand_ratios'] + metal_idx = detailed_dict['metal_idx'] + lig1_atoms = detailed_dict['lig1_atoms'] + lig1_catoms = detailed_dict['lig1_catoms'] + lig2_atoms = detailed_dict['lig2_atoms'] + lig2_catoms = detailed_dict['lig2_catoms'] + lig3_atoms = detailed_dict['lig3_atoms'] + lig3_catoms = detailed_dict['lig3_catoms'] + + if num_unique_ligands == 1: + raise ValueError('Function not defined for homoleptic complexes') + + # two unique ligands: consider cis/trans and fac/mer conversions + elif num_unique_ligands == 2: + if unique_ligand_ratios == 5: + raise ValueError('Function not defined for monoheteroleptic complexes') + + # cis/trans conversion + elif unique_ligand_ratios == 2: + # idx of ligand type 1 to be swapped + # only relevant for cis-->trans conversion, since trans-->cis can be accomplished by switching any two distinct ligands + swap_idx_1 = np.argmin( + np.abs(np.array((tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig2_catoms[1]), + tmc_mol.getAngle(idx0=lig1_catoms[1], idx1=metal_idx, idx2=lig2_catoms[1]), + tmc_mol.getAngle(idx0=lig1_catoms[2], idx1=metal_idx, idx2=lig2_catoms[1]), + tmc_mol.getAngle(idx0=lig1_catoms[3], idx1=metal_idx, idx2=lig2_catoms[1]))) - 180)) + # idx of ligand type 2 to be swapped + # swapping any ligand2 is valid as long as ligand1 is chosen to be orthogonal to it + swap_idx_2 = 0 + + # fac/mer conversion + elif unique_ligand_ratios == 1: + if symmetry == 'mer': + # idx of ligand type 1 to be swapped + swap_idx_1 = np.argmin( + np.abs(np.array((tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig1_catoms[1]), + tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig1_catoms[2]))) - 180)) + 1 + # idx of ligand type 2 to be swapped + swap_idx_2 = np.argmin( + np.abs(np.array((tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig2_catoms[1]), + tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig2_catoms[2]))) - 180)) + 1 + elif symmetry == 'fac': + # idx of ligand type 1 to be swapped + swap_idx_1 = np.argmin( + np.abs(np.array((tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig2_catoms[0]), + tmc_mol.getAngle(idx0=lig1_catoms[1], idx1=metal_idx, idx2=lig2_catoms[0]), + tmc_mol.getAngle(idx0=lig1_catoms[2], idx1=metal_idx, idx2=lig2_catoms[0]))) - 90)) + # idx of ligand type 2 to be swapped + # swapping any ligand2 is valid as long as ligand1 is chosen to be orthogonal to it + swap_idx_2 = 0 + + vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec2 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + atoms_to_move = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] + vec_reflect = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) + vec_reflect = vec_reflect / np.linalg.norm(vec_reflect, 2) + for atom_idx in atoms_to_move: + vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array( + tmc_mol.atoms[metal_idx].coords()) + vec_proj = np.dot(vec_atoms, vec_reflect) * vec_reflect + reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms + tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + + # three unique ligands: consider cis asymmetric (CA)/trans asymmetric (TA), + # double cis symmetric (DCS)/ double trans symmetric (DTS)/equatorial asymmetric (EA), + # and fac asymmetric (FA)/mer asymmetric trans (MAT)/mer asymmetric cis (MAC) conversions + elif num_unique_ligands == 3: + # CA/TA conversion + if unique_ligand_ratios == [4, 1, 4]: + # idx of ligand type 1 to be swapped + # only relevant for CA-->TA conversion, since TA-->CA can be accomplished by switching any two distinct ligands + swap_idx_1 = np.argmin( + np.abs(np.array((tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig3_catoms[0]), + tmc_mol.getAngle(idx0=lig1_catoms[1], idx1=metal_idx, idx2=lig3_catoms[0]), + tmc_mol.getAngle(idx0=lig1_catoms[2], idx1=metal_idx, idx2=lig3_catoms[0]), + tmc_mol.getAngle(idx0=lig1_catoms[3], idx1=metal_idx, idx2=lig3_catoms[0]))) - 180)) + # idx of ligand type 2 to be swapped + swap_idx_2 = 0 + vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec2 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + atoms_to_move = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] + vec_reflect = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) + vec_reflect = vec_reflect / np.linalg.norm(vec_reflect, 2) + for atom_idx in atoms_to_move: + vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array( + tmc_mol.atoms[metal_idx].coords()) + vec_proj = np.dot(vec_atoms, vec_reflect) * vec_reflect + reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms + tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + + # DCS/DTS/EA conversion + elif unique_ligand_ratios == [1, 1, 1]: + # warning: moves L3 to axial position by default, using ligand sorting order from get_symmetry() + if target_symmetry not in ['DCS', 'DTS', 'EA']: + raise ValueError("target_symmetry must be either 'DCS', 'DTS', or 'EA' for TMCs with 3 unique ligands in given stoichiometry") + elif target_symmetry == 'EA' and symmetry == 'double trans symmetric': + # idx of ligand type 1 to be swapped + swap_idx_1 = 0 + # idx of ligand type 2 to be swapped + swap_idx_2 = 0 + vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec2 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + atoms_to_move = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] + vec_reflect = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) + vec_reflect = vec_reflect / np.linalg.norm(vec_reflect, 2) + for atom_idx in atoms_to_move: + vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array( + tmc_mol.atoms[metal_idx].coords()) + vec_proj = np.dot(vec_atoms, vec_reflect) * vec_reflect + reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms + tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + + elif target_symmetry == 'EA' and symmetry == 'double cis symmetric': + # warning: moves L3 to axial position by default, using ligand sorting order from get_symmetry() + # idx of ligand type 2 to be swapped (that which forms 90° angles with both ligands of type 1) + swap_idx_2 = np.argmin(np.array(( + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig2_catoms[0]), + tmc_mol.getAngle(idx0=lig1_catoms[1], idx1=metal_idx, idx2=lig2_catoms[0]))) - 90)), + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig2_catoms[1]), + tmc_mol.getAngle(idx0=lig1_catoms[1], idx1=metal_idx, idx2=lig2_catoms[1]))) - 90))))) + # idx of ligand type 3 to be swapped (that which forms 90° angles with both ligands of type 2) + swap_idx_3 = np.argmin(np.array(( + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig3_catoms[0]), + tmc_mol.getAngle(idx0=lig2_catoms[1], idx1=metal_idx, idx2=lig3_catoms[0]))) - 90)), + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig3_catoms[1]), + tmc_mol.getAngle(idx0=lig2_catoms[1], idx1=metal_idx, idx2=lig3_catoms[1]))) - 90))))) + vec1 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec2 = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + atoms_to_move = lig2_atoms[swap_idx_2] + lig3_atoms[swap_idx_3] + vec_reflect = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) + vec_reflect = vec_reflect / np.linalg.norm(vec_reflect, 2) + for atom_idx in atoms_to_move: + vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec_proj = np.dot(vec_atoms, vec_reflect) * vec_reflect + reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms + tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + + elif target_symmetry == 'DTS' and symmetry == 'equatorial asymmetric': + # figure out which ligands are axial + axial_idx = np.argmin(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig1_catoms[1]), + tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig2_catoms[1]), + tmc_mol.getAngle(idx0=lig3_catoms[0], idx1=metal_idx, idx2=lig3_catoms[1]))) - 180)) + # figure out which ligands are not axial + equatorial_idx = [val for val in range(3) if val != axial_idx] + all_ligand_catoms = [lig1_catoms, lig2_catoms, lig3_catoms] + all_ligand_atoms = [lig1_atoms, lig2_atoms, lig3_atoms] + # idx of equatorial ligand 1 to be swapped. choose arbitarily then choose swap_idx_2 based on this + swap_idx_1 = 0 + # idx of equatorial ligand 2 to be swapped + swap_idx_2 = np.argmin(np.abs(np.array(( + tmc_mol.getAngle(idx0=all_ligand_catoms[equatorial_idx[0]][swap_idx_1], idx1=metal_idx, + idx2=all_ligand_catoms[equatorial_idx[1]][0]), + tmc_mol.getAngle(idx0=all_ligand_catoms[equatorial_idx[0]][swap_idx_1], idx1=metal_idx, + idx2=all_ligand_catoms[equatorial_idx[1]][1]))) - 90)) + + vec1 = np.array(tmc_mol.atoms[all_ligand_catoms[equatorial_idx[0]][swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec2 = np.array(tmc_mol.atoms[all_ligand_catoms[equatorial_idx[1]][swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + atoms_to_move = all_ligand_atoms[equatorial_idx[0]][swap_idx_1] + all_ligand_atoms[equatorial_idx[1]][swap_idx_2] + vec_reflect = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) + vec_reflect = vec_reflect / np.linalg.norm(vec_reflect, 2) + for atom_idx in atoms_to_move: + vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec_proj = np.dot(vec_atoms, vec_reflect) * vec_reflect + reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms + tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + + elif target_symmetry == 'DTS' and symmetry == 'double cis symmetric': + # idx of ligand type 1 to be swapped (that which forms 90° angles with both ligands of type 3) + swap_idx_1 = np.argmin(np.array(( + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig3_catoms[0]), + tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig3_catoms[1]))) - 90)), + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig1_catoms[1], idx1=metal_idx, idx2=lig3_catoms[0]), + tmc_mol.getAngle(idx0=lig1_catoms[1], idx1=metal_idx, idx2=lig3_catoms[1]))) - 90))))) + # idx of ligand type 2 to be swapped (that which forms 90° angles with both ligands of type 1) + swap_idx_2 = np.argmin(np.array(( + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig1_catoms[0]), + tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig1_catoms[1]))) - 90)), + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig2_catoms[1], idx1=metal_idx, idx2=lig1_catoms[0]), + tmc_mol.getAngle(idx0=lig2_catoms[1], idx1=metal_idx, idx2=lig1_catoms[1]))) - 90))))) + # idx of ligand type 3 to be swapped (that which forms 90° angles with both ligands of type 2) + swap_idx_3 = np.argmin(np.array(( + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig3_catoms[0], idx1=metal_idx, idx2=lig2_catoms[0]), + tmc_mol.getAngle(idx0=lig3_catoms[0], idx1=metal_idx, idx2=lig2_catoms[1]))) - 90)), + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig3_catoms[1], idx1=metal_idx, idx2=lig2_catoms[0]), + tmc_mol.getAngle(idx0=lig3_catoms[1], idx1=metal_idx, idx2=lig2_catoms[1]))) - 90))))) + # first reflection + vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec2 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + atoms_to_move_a = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] + vec_reflect_a = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) + vec_reflect_a = vec_reflect_a / np.linalg.norm(vec_reflect_a, 2) + for atom_idx in atoms_to_move_a: + vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec_proj = np.dot(vec_atoms, vec_reflect_a) * vec_reflect_a + reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms + tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + # second reflection + vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec3 = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + atoms_to_move_b = lig1_atoms[swap_idx_1] + lig3_atoms[swap_idx_3] + vec_reflect_b = vec1 / np.linalg.norm(vec1, 2) + vec3 / np.linalg.norm(vec3, 2) + vec_reflect_b = vec_reflect_b / np.linalg.norm(vec_reflect_b, 2) + for atom_idx in atoms_to_move_b: + vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec_proj = np.dot(vec_atoms, vec_reflect_b) * vec_reflect_b + reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms + tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + + elif target_symmetry == 'DCS' and symmetry == 'double trans symmetric': + # these operations can result in unique chiral structures + # idx of ligand type 1 to be swapped + swap_idx_1 = 0 + # idx of ligand type 2 to be swapped + swap_idx_2 = 0 + # idx of ligand type 3 to be swapped + swap_idx_3 = 0 + # first reflection + vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec2 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + atoms_to_move_a = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] + vec_reflect_a = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) + vec_reflect_a = vec_reflect_a / np.linalg.norm(vec_reflect_a, 2) + for atom_idx in atoms_to_move_a: + vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec_proj = np.dot(vec_atoms, vec_reflect_a) * vec_reflect_a + reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms + tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + # second reflection + vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec3 = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + atoms_to_move_b = lig1_atoms[swap_idx_1] + lig3_atoms[swap_idx_3] + vec_reflect_b = vec1 / np.linalg.norm(vec1, 2) + vec3 / np.linalg.norm(vec3, 2) + vec_reflect_b = vec_reflect_b / np.linalg.norm(vec_reflect_b, 2) + for atom_idx in atoms_to_move_b: + vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec_proj = np.dot(vec_atoms, vec_reflect_b) * vec_reflect_b + reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms + tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + + elif target_symmetry == 'DCS' and symmetry == 'equatorial asymmetric': + # these operations can result in unique chiral structures + # figure out which ligands are axial + axial_idx = np.argmin(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig1_catoms[1]), + tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig2_catoms[1]), + tmc_mol.getAngle(idx0=lig3_catoms[0], idx1=metal_idx, idx2=lig3_catoms[1]))) - 180)) + # figure out which ligands are not axial + equatorial_idx = [val for val in range(3) if val != axial_idx] + all_ligand_catoms = [lig1_catoms, lig2_catoms, lig3_catoms] + all_ligand_atoms = [lig1_atoms, lig2_atoms, lig3_atoms] + vec1 = np.array(tmc_mol.atoms[all_ligand_catoms[axial_idx][0]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec2 = np.array(tmc_mol.atoms[all_ligand_catoms[equatorial_idx[0]][0]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + atoms_to_move = all_ligand_atoms[axial_idx][0] + all_ligand_atoms[equatorial_idx[0]][0] + vec_reflect = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) + vec_reflect = vec_reflect / np.linalg.norm(vec_reflect, 2) + for atom_idx in atoms_to_move: + vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec_proj = np.dot(vec_atoms, vec_reflect) * vec_reflect + reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms + tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + + # FA/MAT/MAC conversion + elif unique_ligand_ratios == [3 / 2, 2, 3]: + if target_symmetry not in ['FA', 'MAT', 'MAC']: + raise ValueError("target_symmetry must be either 'FA', 'MAT', or 'MAC' for TMCs with 3 unique ligands in given stoichiometry") + elif target_symmetry == 'FA' and symmetry == 'mer asymmetric trans': + # idx of ligand type 1 to be swapped + swap_idx_1 = np.argmax(np.array(( + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig1_catoms[1]), + tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig1_catoms[2]))) - 90)), + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig1_catoms[1], idx1=metal_idx, idx2=lig1_catoms[0]), + tmc_mol.getAngle(idx0=lig1_catoms[1], idx1=metal_idx, idx2=lig1_catoms[2]))) - 90)), + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig1_catoms[2], idx1=metal_idx, idx2=lig1_catoms[0]), + tmc_mol.getAngle(idx0=lig1_catoms[2], idx1=metal_idx, idx2=lig1_catoms[1]))) - 90))))) + # idx of ligand type 2 to be swapped + swap_idx_2 = 0 + vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec2 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + atoms_to_move = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] + + elif target_symmetry == 'FA' and symmetry == 'mer asymmetric cis': + # idx of ligand type 1 to be swapped + swap_idx_1 = np.argmax(np.array(( + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig1_catoms[1]), + tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig1_catoms[2]))) - 90)), + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig1_catoms[1], idx1=metal_idx, idx2=lig1_catoms[0]), + tmc_mol.getAngle(idx0=lig1_catoms[1], idx1=metal_idx, idx2=lig1_catoms[2]))) - 90)), + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig1_catoms[2], idx1=metal_idx, idx2=lig1_catoms[0]), + tmc_mol.getAngle(idx0=lig1_catoms[2], idx1=metal_idx, idx2=lig1_catoms[1]))) - 90))))) + # idx of ligand type 2 to be swapped + swap_idx_2 = np.argmin(np.array(( + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig1_catoms[0]), + tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig1_catoms[1]), + tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig1_catoms[2]))) - 90)), + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig2_catoms[1], idx1=metal_idx, idx2=lig1_catoms[0]), + tmc_mol.getAngle(idx0=lig2_catoms[1], idx1=metal_idx, idx2=lig1_catoms[1]), + tmc_mol.getAngle(idx0=lig2_catoms[1], idx1=metal_idx, idx2=lig1_catoms[2]))) - 90))))) + vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec2 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + atoms_to_move = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] + + elif target_symmetry == 'MAT' and symmetry == 'fac asymmetric': + # idx of ligand type 1 to be swapped + swap_idx_1 = np.argmin(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig2_catoms[0]), + tmc_mol.getAngle(idx0=lig1_catoms[1], idx1=metal_idx, idx2=lig2_catoms[0]), + tmc_mol.getAngle(idx0=lig1_catoms[2], idx1=metal_idx, idx2=lig2_catoms[0]))) - 180)) + # idx of ligand type 2 to be swapped + swap_idx_2 = 1 + vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec2 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + atoms_to_move = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] + + elif target_symmetry == 'MAT' and symmetry == 'mer asymmetric cis': + # idx of ligand type 2 to be swapped + swap_idx_2 = np.argmax(np.array( + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig1_catoms[0]), + tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig1_catoms[1]), + tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig1_catoms[2]))) - 90)), + np.average(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig2_catoms[1], idx1=metal_idx, idx2=lig1_catoms[0]), + tmc_mol.getAngle(idx0=lig2_catoms[1], idx1=metal_idx, idx2=lig1_catoms[1]), + tmc_mol.getAngle(idx0=lig2_catoms[1], idx1=metal_idx, idx2=lig1_catoms[2]))) - 90)))) + # idx of ligand type 3 to be swapped + swap_idx_3 = 0 + vec1 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec2 = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + atoms_to_move = lig2_atoms[swap_idx_2] + lig3_atoms[swap_idx_3] + + elif target_symmetry == 'MAC' and symmetry == 'fac asymmetric': + # these operations can result in unique chiral structures + # idx of ligand type 1 to be swapped + swap_idx_1 = np.argmin(np.abs(np.array(( + tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig3_catoms[0]), + tmc_mol.getAngle(idx0=lig1_catoms[1], idx1=metal_idx, idx2=lig3_catoms[0]), + tmc_mol.getAngle(idx0=lig1_catoms[2], idx1=metal_idx, idx2=lig3_catoms[0]))) - 90)) + # idx of ligand type 3 to be swapped + swap_idx_3 = 0 + vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec2 = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + atoms_to_move = lig1_atoms[swap_idx_1] + lig3_atoms[swap_idx_3] + + elif target_symmetry == 'MAC' and symmetry == 'mer asymmetric trans': + # idx of ligand type 2 to be swapped + swap_idx_2 = 0 + # idx of ligand type 3 to be swapped + swap_idx_3 = 0 + vec1 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec2 = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + atoms_to_move = lig2_atoms[swap_idx_2] + lig3_atoms[swap_idx_3] + + vec_reflect = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) + vec_reflect = vec_reflect / np.linalg.norm(vec_reflect, 2) + for atom_idx in atoms_to_move: + vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + vec_proj = np.dot(vec_atoms, vec_reflect) * vec_reflect + reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms + tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + + return tmc_mol def get_features(self, lac=True, force_generate=False, eq_sym=False, use_dist=False, NumB=False, Gval=False, size_normalize=False, From 703a06bf5a646267b005c08fdafbf6cf2c9e6016 Mon Sep 17 00:00:00 2001 From: jwtoney Date: Wed, 11 Dec 2024 16:11:37 -0500 Subject: [PATCH 2/3] added helper functions, improved documentation --- molSimplify/Classes/mol3D.py | 229 +++++++++++++++++------------------ 1 file changed, 110 insertions(+), 119 deletions(-) diff --git a/molSimplify/Classes/mol3D.py b/molSimplify/Classes/mol3D.py index 4d69f1da..2259d262 100644 --- a/molSimplify/Classes/mol3D.py +++ b/molSimplify/Classes/mol3D.py @@ -6372,18 +6372,13 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal # idx of ligand type 2 to be swapped # swapping any ligand2 is valid as long as ligand1 is chosen to be orthogonal to it swap_idx_2 = 0 - - vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec2 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) + lig2_catom_coords = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) + # define all atoms to be moved, i.e., all atoms in both ligands that will be moved atoms_to_move = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] - vec_reflect = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) - vec_reflect = vec_reflect / np.linalg.norm(vec_reflect, 2) - for atom_idx in atoms_to_move: - vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array( - tmc_mol.atoms[metal_idx].coords()) - vec_proj = np.dot(vec_atoms, vec_reflect) * vec_reflect - reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms - tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), + lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, + atoms_to_move=atoms_to_move) # three unique ligands: consider cis asymmetric (CA)/trans asymmetric (TA), # double cis symmetric (DCS)/ double trans symmetric (DTS)/equatorial asymmetric (EA), @@ -6400,17 +6395,12 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.getAngle(idx0=lig1_catoms[3], idx1=metal_idx, idx2=lig3_catoms[0]))) - 180)) # idx of ligand type 2 to be swapped swap_idx_2 = 0 - vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec2 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) + lig2_catom_coords = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) atoms_to_move = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] - vec_reflect = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) - vec_reflect = vec_reflect / np.linalg.norm(vec_reflect, 2) - for atom_idx in atoms_to_move: - vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array( - tmc_mol.atoms[metal_idx].coords()) - vec_proj = np.dot(vec_atoms, vec_reflect) * vec_reflect - reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms - tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), + lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, + atoms_to_move=atoms_to_move) # DCS/DTS/EA conversion elif unique_ligand_ratios == [1, 1, 1]: @@ -6422,17 +6412,12 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal swap_idx_1 = 0 # idx of ligand type 2 to be swapped swap_idx_2 = 0 - vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec2 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) + lig2_catom_coords = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) atoms_to_move = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] - vec_reflect = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) - vec_reflect = vec_reflect / np.linalg.norm(vec_reflect, 2) - for atom_idx in atoms_to_move: - vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array( - tmc_mol.atoms[metal_idx].coords()) - vec_proj = np.dot(vec_atoms, vec_reflect) * vec_reflect - reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms - tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), + lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, + atoms_to_move=atoms_to_move) elif target_symmetry == 'EA' and symmetry == 'double cis symmetric': # warning: moves L3 to axial position by default, using ligand sorting order from get_symmetry() @@ -6452,16 +6437,12 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal np.average(np.abs(np.array(( tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig3_catoms[1]), tmc_mol.getAngle(idx0=lig2_catoms[1], idx1=metal_idx, idx2=lig3_catoms[1]))) - 90))))) - vec1 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec2 = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + lig1_catom_coords = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) + lig2_catom_coords = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) atoms_to_move = lig2_atoms[swap_idx_2] + lig3_atoms[swap_idx_3] - vec_reflect = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) - vec_reflect = vec_reflect / np.linalg.norm(vec_reflect, 2) - for atom_idx in atoms_to_move: - vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec_proj = np.dot(vec_atoms, vec_reflect) * vec_reflect - reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms - tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), + lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, + atoms_to_move=atoms_to_move) elif target_symmetry == 'DTS' and symmetry == 'equatorial asymmetric': # figure out which ligands are axial @@ -6482,16 +6463,12 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.getAngle(idx0=all_ligand_catoms[equatorial_idx[0]][swap_idx_1], idx1=metal_idx, idx2=all_ligand_catoms[equatorial_idx[1]][1]))) - 90)) - vec1 = np.array(tmc_mol.atoms[all_ligand_catoms[equatorial_idx[0]][swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec2 = np.array(tmc_mol.atoms[all_ligand_catoms[equatorial_idx[1]][swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + lig1_catom_coords = np.array(tmc_mol.atoms[all_ligand_catoms[equatorial_idx[0]][swap_idx_1]].coords()) + lig2_catom_coords = np.array(tmc_mol.atoms[all_ligand_catoms[equatorial_idx[1]][swap_idx_2]].coords()) atoms_to_move = all_ligand_atoms[equatorial_idx[0]][swap_idx_1] + all_ligand_atoms[equatorial_idx[1]][swap_idx_2] - vec_reflect = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) - vec_reflect = vec_reflect / np.linalg.norm(vec_reflect, 2) - for atom_idx in atoms_to_move: - vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec_proj = np.dot(vec_atoms, vec_reflect) * vec_reflect - reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms - tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), + lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, + atoms_to_move=atoms_to_move) elif target_symmetry == 'DTS' and symmetry == 'double cis symmetric': # idx of ligand type 1 to be swapped (that which forms 90° angles with both ligands of type 3) @@ -6519,27 +6496,19 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.getAngle(idx0=lig3_catoms[1], idx1=metal_idx, idx2=lig2_catoms[0]), tmc_mol.getAngle(idx0=lig3_catoms[1], idx1=metal_idx, idx2=lig2_catoms[1]))) - 90))))) # first reflection - vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec2 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - atoms_to_move_a = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] - vec_reflect_a = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) - vec_reflect_a = vec_reflect_a / np.linalg.norm(vec_reflect_a, 2) - for atom_idx in atoms_to_move_a: - vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec_proj = np.dot(vec_atoms, vec_reflect_a) * vec_reflect_a - reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms - tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) + lig2_catom_coords = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) + atoms_to_move = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] + tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), + lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, + atoms_to_move=atoms_to_move) # second reflection - vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec3 = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - atoms_to_move_b = lig1_atoms[swap_idx_1] + lig3_atoms[swap_idx_3] - vec_reflect_b = vec1 / np.linalg.norm(vec1, 2) + vec3 / np.linalg.norm(vec3, 2) - vec_reflect_b = vec_reflect_b / np.linalg.norm(vec_reflect_b, 2) - for atom_idx in atoms_to_move_b: - vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec_proj = np.dot(vec_atoms, vec_reflect_b) * vec_reflect_b - reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms - tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) + lig2_catom_coords = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) + atoms_to_move = lig1_atoms[swap_idx_1] + lig3_atoms[swap_idx_3] + tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), + lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, + atoms_to_move=atoms_to_move) elif target_symmetry == 'DCS' and symmetry == 'double trans symmetric': # these operations can result in unique chiral structures @@ -6550,27 +6519,33 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal # idx of ligand type 3 to be swapped swap_idx_3 = 0 # first reflection - vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec2 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - atoms_to_move_a = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] - vec_reflect_a = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) - vec_reflect_a = vec_reflect_a / np.linalg.norm(vec_reflect_a, 2) - for atom_idx in atoms_to_move_a: - vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec_proj = np.dot(vec_atoms, vec_reflect_a) * vec_reflect_a - reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms - tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) + lig2_catom_coords = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) + atoms_to_move = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] + tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), + lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, + atoms_to_move=atoms_to_move) + # vec_reflect_a = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) + # vec_reflect_a = vec_reflect_a / np.linalg.norm(vec_reflect_a, 2) + # for atom_idx in atoms_to_move_a: + # vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + # vec_proj = np.dot(vec_atoms, vec_reflect_a) * vec_reflect_a + # reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms + # tmc_mol.atoms[atom_idx].setcoords(reflected_coords) # second reflection - vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec3 = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - atoms_to_move_b = lig1_atoms[swap_idx_1] + lig3_atoms[swap_idx_3] - vec_reflect_b = vec1 / np.linalg.norm(vec1, 2) + vec3 / np.linalg.norm(vec3, 2) - vec_reflect_b = vec_reflect_b / np.linalg.norm(vec_reflect_b, 2) - for atom_idx in atoms_to_move_b: - vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec_proj = np.dot(vec_atoms, vec_reflect_b) * vec_reflect_b - reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms - tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) + lig2_catom_coords = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) + atoms_to_move = lig1_atoms[swap_idx_1] + lig3_atoms[swap_idx_3] + tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), + lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, + atoms_to_move=atoms_to_move) + # vec_reflect_b = vec1 / np.linalg.norm(vec1, 2) + vec3 / np.linalg.norm(vec3, 2) + # vec_reflect_b = vec_reflect_b / np.linalg.norm(vec_reflect_b, 2) + # for atom_idx in atoms_to_move_b: + # vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + # vec_proj = np.dot(vec_atoms, vec_reflect_b) * vec_reflect_b + # reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms + # tmc_mol.atoms[atom_idx].setcoords(reflected_coords) elif target_symmetry == 'DCS' and symmetry == 'equatorial asymmetric': # these operations can result in unique chiral structures @@ -6583,16 +6558,12 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal equatorial_idx = [val for val in range(3) if val != axial_idx] all_ligand_catoms = [lig1_catoms, lig2_catoms, lig3_catoms] all_ligand_atoms = [lig1_atoms, lig2_atoms, lig3_atoms] - vec1 = np.array(tmc_mol.atoms[all_ligand_catoms[axial_idx][0]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec2 = np.array(tmc_mol.atoms[all_ligand_catoms[equatorial_idx[0]][0]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + lig1_catom_coords = np.array(tmc_mol.atoms[all_ligand_catoms[axial_idx][0]].coords()) + lig2_catom_coords = np.array(tmc_mol.atoms[all_ligand_catoms[equatorial_idx[0]][0]].coords()) atoms_to_move = all_ligand_atoms[axial_idx][0] + all_ligand_atoms[equatorial_idx[0]][0] - vec_reflect = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) - vec_reflect = vec_reflect / np.linalg.norm(vec_reflect, 2) - for atom_idx in atoms_to_move: - vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec_proj = np.dot(vec_atoms, vec_reflect) * vec_reflect - reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms - tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), + lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, + atoms_to_move=atoms_to_move) # FA/MAT/MAC conversion elif unique_ligand_ratios == [3 / 2, 2, 3]: @@ -6612,8 +6583,8 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.getAngle(idx0=lig1_catoms[2], idx1=metal_idx, idx2=lig1_catoms[1]))) - 90))))) # idx of ligand type 2 to be swapped swap_idx_2 = 0 - vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec2 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) + lig2_catom_coords = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) atoms_to_move = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] elif target_symmetry == 'FA' and symmetry == 'mer asymmetric cis': @@ -6638,8 +6609,8 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.getAngle(idx0=lig2_catoms[1], idx1=metal_idx, idx2=lig1_catoms[0]), tmc_mol.getAngle(idx0=lig2_catoms[1], idx1=metal_idx, idx2=lig1_catoms[1]), tmc_mol.getAngle(idx0=lig2_catoms[1], idx1=metal_idx, idx2=lig1_catoms[2]))) - 90))))) - vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec2 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) + lig2_catom_coords = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) atoms_to_move = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] elif target_symmetry == 'MAT' and symmetry == 'fac asymmetric': @@ -6650,8 +6621,8 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.getAngle(idx0=lig1_catoms[2], idx1=metal_idx, idx2=lig2_catoms[0]))) - 180)) # idx of ligand type 2 to be swapped swap_idx_2 = 1 - vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec2 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) + lig2_catom_coords = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) atoms_to_move = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] elif target_symmetry == 'MAT' and symmetry == 'mer asymmetric cis': @@ -6667,8 +6638,8 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.getAngle(idx0=lig2_catoms[1], idx1=metal_idx, idx2=lig1_catoms[2]))) - 90)))) # idx of ligand type 3 to be swapped swap_idx_3 = 0 - vec1 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec2 = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + lig1_catom_coords = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) + lig2_catom_coords = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) atoms_to_move = lig2_atoms[swap_idx_2] + lig3_atoms[swap_idx_3] elif target_symmetry == 'MAC' and symmetry == 'fac asymmetric': @@ -6680,8 +6651,8 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.getAngle(idx0=lig1_catoms[2], idx1=metal_idx, idx2=lig3_catoms[0]))) - 90)) # idx of ligand type 3 to be swapped swap_idx_3 = 0 - vec1 = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec2 = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) + lig2_catom_coords = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) atoms_to_move = lig1_atoms[swap_idx_1] + lig3_atoms[swap_idx_3] elif target_symmetry == 'MAC' and symmetry == 'mer asymmetric trans': @@ -6689,18 +6660,38 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal swap_idx_2 = 0 # idx of ligand type 3 to be swapped swap_idx_3 = 0 - vec1 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec2 = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + # vec1 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + # vec2 = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + lig1_catom_coords = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) + lig2_catom_coords = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) atoms_to_move = lig2_atoms[swap_idx_2] + lig3_atoms[swap_idx_3] - - vec_reflect = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) - vec_reflect = vec_reflect / np.linalg.norm(vec_reflect, 2) - for atom_idx in atoms_to_move: - vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - vec_proj = np.dot(vec_atoms, vec_reflect) * vec_reflect - reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms - tmc_mol.atoms[atom_idx].setcoords(reflected_coords) - + tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), + lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, + atoms_to_move=atoms_to_move) + + # vec_reflect = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) + # vec_reflect = vec_reflect / np.linalg.norm(vec_reflect, 2) + # for atom_idx in atoms_to_move: + # vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) + # vec_proj = np.dot(vec_atoms, vec_reflect) * vec_reflect + # reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms + # tmc_mol.atoms[atom_idx].setcoords(reflected_coords) + return tmc_mol + def reflect_coords(tmc_mol, metal_coords, lig1_catom_coords, lig2_catom_coords, atoms_to_move): + # define vector from metal to first coordinating atom of first ligand to be swapped + vec1 = lig1_catom_coords - metal_coords + # define vector from metal to first coordinating atom of second ligand to be swapped + vec2 = lig2_catom_coords - metal_coords + # define reflection vector between the two previous vectors and normalize + vec_reflect = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) + vec_reflect = vec_reflect / np.linalg.norm(vec_reflect, 2) + # define all atoms to be moved, i.e., all atoms in both ligands that will be moved + # move all atoms by flipping coordinates across normalized reflection vector + for atom_idx in atoms_to_move: + vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - metal_coords + vec_proj = np.dot(vec_atoms, vec_reflect) * vec_reflect + reflected_coords = metal_coords + 2 * vec_proj - vec_atoms + tmc_mol.atoms[atom_idx].setcoords(reflected_coords) return tmc_mol def get_features(self, lac=True, force_generate=False, eq_sym=False, From a44f4cea6b887351be94808075c275e56fbe2ddb Mon Sep 17 00:00:00 2001 From: jwtoney Date: Wed, 11 Dec 2024 16:21:23 -0500 Subject: [PATCH 3/3] added helper functions, improved documentation --- molSimplify/Classes/mol3D.py | 95 ++++++++++++++++-------------------- 1 file changed, 41 insertions(+), 54 deletions(-) diff --git a/molSimplify/Classes/mol3D.py b/molSimplify/Classes/mol3D.py index 2259d262..fcc49055 100644 --- a/molSimplify/Classes/mol3D.py +++ b/molSimplify/Classes/mol3D.py @@ -6320,7 +6320,6 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal max_allowed_dev=max_allowed_dev, details=True) symmetry = symmetry_dict['symmetry'] - num_unique_ligands = detailed_dict['num_unique_ligands'] unique_ligand_ratios = detailed_dict['unique_ligand_ratios'] metal_idx = detailed_dict['metal_idx'] @@ -6333,12 +6332,10 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal if num_unique_ligands == 1: raise ValueError('Function not defined for homoleptic complexes') - # two unique ligands: consider cis/trans and fac/mer conversions elif num_unique_ligands == 2: if unique_ligand_ratios == 5: raise ValueError('Function not defined for monoheteroleptic complexes') - # cis/trans conversion elif unique_ligand_ratios == 2: # idx of ligand type 1 to be swapped @@ -6348,10 +6345,8 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.getAngle(idx0=lig1_catoms[1], idx1=metal_idx, idx2=lig2_catoms[1]), tmc_mol.getAngle(idx0=lig1_catoms[2], idx1=metal_idx, idx2=lig2_catoms[1]), tmc_mol.getAngle(idx0=lig1_catoms[3], idx1=metal_idx, idx2=lig2_catoms[1]))) - 180)) - # idx of ligand type 2 to be swapped - # swapping any ligand2 is valid as long as ligand1 is chosen to be orthogonal to it + # idx of ligand type 2 to be swapped (any ligand2 is valid as long as ligand1 is orthogonal) swap_idx_2 = 0 - # fac/mer conversion elif unique_ligand_ratios == 1: if symmetry == 'mer': @@ -6369,8 +6364,7 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal np.abs(np.array((tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig2_catoms[0]), tmc_mol.getAngle(idx0=lig1_catoms[1], idx1=metal_idx, idx2=lig2_catoms[0]), tmc_mol.getAngle(idx0=lig1_catoms[2], idx1=metal_idx, idx2=lig2_catoms[0]))) - 90)) - # idx of ligand type 2 to be swapped - # swapping any ligand2 is valid as long as ligand1 is chosen to be orthogonal to it + # idx of ligand type 2 to be swapped (ligand2 is valid as long as ligand1 is orthogonal) swap_idx_2 = 0 lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) lig2_catom_coords = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) @@ -6379,7 +6373,6 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, atoms_to_move=atoms_to_move) - # three unique ligands: consider cis asymmetric (CA)/trans asymmetric (TA), # double cis symmetric (DCS)/ double trans symmetric (DTS)/equatorial asymmetric (EA), # and fac asymmetric (FA)/mer asymmetric trans (MAT)/mer asymmetric cis (MAC) conversions @@ -6401,16 +6394,15 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, atoms_to_move=atoms_to_move) - # DCS/DTS/EA conversion elif unique_ligand_ratios == [1, 1, 1]: # warning: moves L3 to axial position by default, using ligand sorting order from get_symmetry() if target_symmetry not in ['DCS', 'DTS', 'EA']: raise ValueError("target_symmetry must be either 'DCS', 'DTS', or 'EA' for TMCs with 3 unique ligands in given stoichiometry") + # DTS to EA elif target_symmetry == 'EA' and symmetry == 'double trans symmetric': - # idx of ligand type 1 to be swapped + # idx of ligand types 1 and 2 to be swapped swap_idx_1 = 0 - # idx of ligand type 2 to be swapped swap_idx_2 = 0 lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) lig2_catom_coords = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) @@ -6418,7 +6410,7 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, atoms_to_move=atoms_to_move) - + # DCS to EA elif target_symmetry == 'EA' and symmetry == 'double cis symmetric': # warning: moves L3 to axial position by default, using ligand sorting order from get_symmetry() # idx of ligand type 2 to be swapped (that which forms 90° angles with both ligands of type 1) @@ -6443,14 +6435,14 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, atoms_to_move=atoms_to_move) - + # EA to DTS elif target_symmetry == 'DTS' and symmetry == 'equatorial asymmetric': # figure out which ligands are axial axial_idx = np.argmin(np.abs(np.array(( tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig1_catoms[1]), tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig2_catoms[1]), tmc_mol.getAngle(idx0=lig3_catoms[0], idx1=metal_idx, idx2=lig3_catoms[1]))) - 180)) - # figure out which ligands are not axial + # figure out which ligands are not axial (i.e., equatorial) equatorial_idx = [val for val in range(3) if val != axial_idx] all_ligand_catoms = [lig1_catoms, lig2_catoms, lig3_catoms] all_ligand_atoms = [lig1_atoms, lig2_atoms, lig3_atoms] @@ -6469,7 +6461,7 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, atoms_to_move=atoms_to_move) - + # DCS to DTS elif target_symmetry == 'DTS' and symmetry == 'double cis symmetric': # idx of ligand type 1 to be swapped (that which forms 90° angles with both ligands of type 3) swap_idx_1 = np.argmin(np.array(( @@ -6509,14 +6501,12 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, atoms_to_move=atoms_to_move) - + # DTS to DCS elif target_symmetry == 'DCS' and symmetry == 'double trans symmetric': # these operations can result in unique chiral structures - # idx of ligand type 1 to be swapped + # idx of ligand types 1, 2, and 3 to be swapped swap_idx_1 = 0 - # idx of ligand type 2 to be swapped swap_idx_2 = 0 - # idx of ligand type 3 to be swapped swap_idx_3 = 0 # first reflection lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) @@ -6525,13 +6515,6 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, atoms_to_move=atoms_to_move) - # vec_reflect_a = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) - # vec_reflect_a = vec_reflect_a / np.linalg.norm(vec_reflect_a, 2) - # for atom_idx in atoms_to_move_a: - # vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - # vec_proj = np.dot(vec_atoms, vec_reflect_a) * vec_reflect_a - # reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms - # tmc_mol.atoms[atom_idx].setcoords(reflected_coords) # second reflection lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) lig2_catom_coords = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) @@ -6539,14 +6522,7 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, atoms_to_move=atoms_to_move) - # vec_reflect_b = vec1 / np.linalg.norm(vec1, 2) + vec3 / np.linalg.norm(vec3, 2) - # vec_reflect_b = vec_reflect_b / np.linalg.norm(vec_reflect_b, 2) - # for atom_idx in atoms_to_move_b: - # vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - # vec_proj = np.dot(vec_atoms, vec_reflect_b) * vec_reflect_b - # reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms - # tmc_mol.atoms[atom_idx].setcoords(reflected_coords) - + # EA to DCS elif target_symmetry == 'DCS' and symmetry == 'equatorial asymmetric': # these operations can result in unique chiral structures # figure out which ligands are axial @@ -6554,7 +6530,7 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.getAngle(idx0=lig1_catoms[0], idx1=metal_idx, idx2=lig1_catoms[1]), tmc_mol.getAngle(idx0=lig2_catoms[0], idx1=metal_idx, idx2=lig2_catoms[1]), tmc_mol.getAngle(idx0=lig3_catoms[0], idx1=metal_idx, idx2=lig3_catoms[1]))) - 180)) - # figure out which ligands are not axial + # figure out which ligands are not axial (i.e., equatorial) equatorial_idx = [val for val in range(3) if val != axial_idx] all_ligand_catoms = [lig1_catoms, lig2_catoms, lig3_catoms] all_ligand_atoms = [lig1_atoms, lig2_atoms, lig3_atoms] @@ -6564,11 +6540,11 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, atoms_to_move=atoms_to_move) - # FA/MAT/MAC conversion elif unique_ligand_ratios == [3 / 2, 2, 3]: if target_symmetry not in ['FA', 'MAT', 'MAC']: raise ValueError("target_symmetry must be either 'FA', 'MAT', or 'MAC' for TMCs with 3 unique ligands in given stoichiometry") + # MAT to FA elif target_symmetry == 'FA' and symmetry == 'mer asymmetric trans': # idx of ligand type 1 to be swapped swap_idx_1 = np.argmax(np.array(( @@ -6586,7 +6562,7 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) lig2_catom_coords = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) atoms_to_move = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] - + # MAC to FA elif target_symmetry == 'FA' and symmetry == 'mer asymmetric cis': # idx of ligand type 1 to be swapped swap_idx_1 = np.argmax(np.array(( @@ -6612,7 +6588,7 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) lig2_catom_coords = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) atoms_to_move = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] - + # FA to MAT elif target_symmetry == 'MAT' and symmetry == 'fac asymmetric': # idx of ligand type 1 to be swapped swap_idx_1 = np.argmin(np.abs(np.array(( @@ -6624,7 +6600,7 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) lig2_catom_coords = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) atoms_to_move = lig1_atoms[swap_idx_1] + lig2_atoms[swap_idx_2] - + # MAC to MAT elif target_symmetry == 'MAT' and symmetry == 'mer asymmetric cis': # idx of ligand type 2 to be swapped swap_idx_2 = np.argmax(np.array( @@ -6641,7 +6617,7 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal lig1_catom_coords = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) lig2_catom_coords = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) atoms_to_move = lig2_atoms[swap_idx_2] + lig3_atoms[swap_idx_3] - + # FA to MAC elif target_symmetry == 'MAC' and symmetry == 'fac asymmetric': # these operations can result in unique chiral structures # idx of ligand type 1 to be swapped @@ -6654,30 +6630,41 @@ def flip_symmetry(tmc_mol, verbose=True, max_allowed_dev=30, target_symmetry=Fal lig1_catom_coords = np.array(tmc_mol.atoms[lig1_catoms[swap_idx_1]].coords()) lig2_catom_coords = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) atoms_to_move = lig1_atoms[swap_idx_1] + lig3_atoms[swap_idx_3] - + # MAT to MAC elif target_symmetry == 'MAC' and symmetry == 'mer asymmetric trans': - # idx of ligand type 2 to be swapped + # idx of ligand types 2 and 3 to be swapped swap_idx_2 = 0 - # idx of ligand type 3 to be swapped swap_idx_3 = 0 - # vec1 = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - # vec2 = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) lig1_catom_coords = np.array(tmc_mol.atoms[lig2_catoms[swap_idx_2]].coords()) lig2_catom_coords = np.array(tmc_mol.atoms[lig3_catoms[swap_idx_3]].coords()) atoms_to_move = lig2_atoms[swap_idx_2] + lig3_atoms[swap_idx_3] tmc_mol.reflect_coords(metal_coords=np.array(tmc_mol.atoms[metal_idx].coords()), lig1_catom_coords=lig1_catom_coords, lig2_catom_coords=lig2_catom_coords, atoms_to_move=atoms_to_move) - - # vec_reflect = vec1 / np.linalg.norm(vec1, 2) + vec2 / np.linalg.norm(vec2, 2) - # vec_reflect = vec_reflect / np.linalg.norm(vec_reflect, 2) - # for atom_idx in atoms_to_move: - # vec_atoms = np.array(tmc_mol.atoms[atom_idx].coords()) - np.array(tmc_mol.atoms[metal_idx].coords()) - # vec_proj = np.dot(vec_atoms, vec_reflect) * vec_reflect - # reflected_coords = np.array(tmc_mol.atoms[metal_idx].coords()) + 2 * vec_proj - vec_atoms - # tmc_mol.atoms[atom_idx].setcoords(reflected_coords) return tmc_mol def reflect_coords(tmc_mol, metal_coords, lig1_catom_coords, lig2_catom_coords, atoms_to_move): + """ + Helper function for flip symmetry to calculate vectors, projections, and update coordinates + + Parameters + ---------- + tmc_mol: mol3D + mol3D instance of TMC + metal_coords: np.array + Coordinates of metal in TMC + lig1_catom_coords: np.array + Coordinates of coordinating atom of first ligand to be flipped + lig2_catom_coords: np.array + Coordinates of coordinating atom of second ligand to be flipped + atoms_to_move: list + List of atom indices to be moved + + Returns + ------- + tmc_mol: mol3D + returns self, a mol3D object with flipped symmetry + """ + # define vector from metal to first coordinating atom of first ligand to be swapped vec1 = lig1_catom_coords - metal_coords # define vector from metal to first coordinating atom of second ligand to be swapped