From f5dc869c4ecccafa2bde637e1f694058d475abc9 Mon Sep 17 00:00:00 2001 From: aarongarrison Date: Wed, 18 Sep 2024 17:33:22 -0400 Subject: [PATCH 1/3] Allowing for easy access of the aligned structures during the align_rmsd program. --- molSimplify/Scripts/rmsd.py | 54 ++++++++++++++++++++++++++++--------- 1 file changed, 42 insertions(+), 12 deletions(-) diff --git a/molSimplify/Scripts/rmsd.py b/molSimplify/Scripts/rmsd.py index fc67b83f..e9f691b3 100644 --- a/molSimplify/Scripts/rmsd.py +++ b/molSimplify/Scripts/rmsd.py @@ -293,7 +293,7 @@ def reorder_hungarian(p_atoms, q_atoms, p_coord, q_coord): ---------- p_atoms : np.array (N,1) matrix, where N is points holding the atoms' names - p_atoms : np.array + q_atoms : np.array (N,1) matrix, where N is points holding the atoms' names p_coord : np.array (N,D) matrix, where N is points and D is dimension @@ -450,7 +450,7 @@ def rmsd_reorder_rotate(p_atoms, q_atoms, p_coord, q_coord, def reorder_rotate(p_atoms, q_atoms, p_coord, q_coord, rotation="kabsch", reorder="hungarian", - translate=True): + translate=True, return_reorder=False): """Reorders atoms pairwise and rotates structures onto one another. Parameters @@ -471,6 +471,9 @@ def reorder_rotate(p_atoms, q_atoms, p_coord, q_coord, Whether or not the molecules should be translated such that their centroid is at the origin. Default is True. + return_reorder : bool, optional + Whether or not the reordering is returned. + Default is False. Returns ------- @@ -516,12 +519,15 @@ def reorder_rotate(p_atoms, q_atoms, p_coord, q_coord, q_review = reorder_method(p_atoms, q_atoms, p_coord, q_coord) q_coord = q_coord[q_review] # q_atoms = q_atoms[q_review] - # print("q_review", q_review) + #print("q_review", q_review) if rotation_method is not None: q_coord = rotation_method(q_coord, p_coord) - return q_coord + if return_reorder: + return q_coord, q_review + else: + return q_coord def rigorous_rmsd(mol_p, mol_q, rotation: str = "kabsch", @@ -679,6 +685,21 @@ def align_rmsd_rotate(mol_p, mol_q, rotation: str = "kabsch", molq_atoms = mol_q.symvect() molq_coords = rotate_onto_principal_axes(mol_q) + if verbose: + print('Molecule 1 on principal axes:') + print(mol_p.natoms, '\n') + for i, atom in enumerate(mol_p.getAtoms()): + xyz = molp_coords[i] + ss = "%s \t%f\t%f\t%f" % (atom.sym, xyz[0], xyz[1], xyz[2]) + print(ss) + + print('Molecule 2 on principal axes:') + print(mol_q.natoms, '\n') + for i, atom in enumerate(mol_q.getAtoms()): + xyz = molq_coords[i] + ss = "%s \t%f\t%f\t%f" % (atom.sym, xyz[0], xyz[1], xyz[2]) + print(ss) + #allow 180 degree rotations about each axis to catch if aligned to different sides best_rmsd = np.inf x_rot = np.array([[1, 0, 0], [0, -1, 0], [0, 0, -1]]) #180 about x @@ -686,18 +707,19 @@ def align_rmsd_rotate(mol_p, mol_q, rotation: str = "kabsch", z_rot = np.array([[-1, 0, 0], [0, -1, 0], [0, 0, 1]]) #180 about z transformations = [ np.eye(3), #no change - x_rot, y_rot, z_rot, - x_rot@y_rot, x_rot@z_rot, y_rot@z_rot, - x_rot@y_rot@z_rot - ] + x_rot, y_rot, z_rot] for transformation in transformations: transformed_molq_coords = np.vstack([transformation @ molq_coords[i, :] for i in range(len(molq_coords))]) #Iterate for the specified number of iterations for i in range(iterations): - transformed_molq_coords = reorder_rotate(molp_atoms, molq_atoms, molp_coords, transformed_molq_coords, - rotation=rotation, reorder=reorder, translate=True) + if verbose: + return_reorder = True + else: + return_reorder = False + transformed_molq_coords, q_review = reorder_rotate(molp_atoms, molq_atoms, molp_coords, transformed_molq_coords, + rotation=rotation, reorder=reorder, translate=True, return_reorder=return_reorder) if i == iterations-1: #for the final iteration, compute the RMSD and compare @@ -705,6 +727,16 @@ def align_rmsd_rotate(mol_p, mol_q, rotation: str = "kabsch", if result_rmsd < best_rmsd: best_rmsd = result_rmsd + if verbose: + print('Aligned molecule 2:') + print(f'The transformation is {transformation}.') + print(f'The RMSD is {result_rmsd}.') + print(mol_q.natoms, '\n') + symbs = molq_atoms[q_review] + for j, atom in enumerate(mol_q.getAtoms()): + xyz = transformed_molq_coords[j] + ss = "%s \t%f\t%f\t%f" % (symbs[j], xyz[0], xyz[1], xyz[2]) + print(ss) if verbose: cutoff = 10 #How close moments have to be for a warning @@ -719,8 +751,6 @@ def align_rmsd_rotate(mol_p, mol_q, rotation: str = "kabsch", print(pmom2) print('This may lead to improper orientation in some cases.') - #Can address improper rotations by also allowing 90 degree rotations, combinations (including 180 rotations as well) - return best_rmsd def test_case(): From d814cee3c06c2089544e37d13682d59c327cb3e6 Mon Sep 17 00:00:00 2001 From: aarongarrison Date: Wed, 18 Sep 2024 17:57:35 -0400 Subject: [PATCH 2/3] Bugfix for verbose option. --- molSimplify/Scripts/rmsd.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/molSimplify/Scripts/rmsd.py b/molSimplify/Scripts/rmsd.py index e9f691b3..fc93273b 100644 --- a/molSimplify/Scripts/rmsd.py +++ b/molSimplify/Scripts/rmsd.py @@ -715,11 +715,12 @@ def align_rmsd_rotate(mol_p, mol_q, rotation: str = "kabsch", #Iterate for the specified number of iterations for i in range(iterations): if verbose: - return_reorder = True + transformed_molq_coords, q_review = reorder_rotate(molp_atoms, molq_atoms, molp_coords, transformed_molq_coords, + rotation=rotation, reorder=reorder, translate=True, return_reorder=True) else: - return_reorder = False - transformed_molq_coords, q_review = reorder_rotate(molp_atoms, molq_atoms, molp_coords, transformed_molq_coords, - rotation=rotation, reorder=reorder, translate=True, return_reorder=return_reorder) + transformed_molq_coords = reorder_rotate(molp_atoms, molq_atoms, molp_coords, transformed_molq_coords, + rotation=rotation, reorder=reorder, translate=True, return_reorder=False) + if i == iterations-1: #for the final iteration, compute the RMSD and compare From 60c89ca75aafe8155ff01b7eb8d191f0c49f45b3 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 18 Sep 2024 22:05:40 +0000 Subject: [PATCH 3/3] [pre-commit.ci] auto fixes from pre-commit hooks --- molSimplify/Scripts/rmsd.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/molSimplify/Scripts/rmsd.py b/molSimplify/Scripts/rmsd.py index fc93273b..4811cf32 100644 --- a/molSimplify/Scripts/rmsd.py +++ b/molSimplify/Scripts/rmsd.py @@ -720,7 +720,7 @@ def align_rmsd_rotate(mol_p, mol_q, rotation: str = "kabsch", else: transformed_molq_coords = reorder_rotate(molp_atoms, molq_atoms, molp_coords, transformed_molq_coords, rotation=rotation, reorder=reorder, translate=True, return_reorder=False) - + if i == iterations-1: #for the final iteration, compute the RMSD and compare