lib.structure.superimpose

1 ############################################################################### 2 # # 3 # Copyright (C) 2011-2013 Edward d'Auvergne # 4 # # 5 # This file is part of the program relax (http://www.nmr-relax.com). # 6 # # 7 # This program is free software: you can redistribute it and/or modify # 8 # it under the terms of the GNU General Public License as published by # 9 # the Free Software Foundation, either version 3 of the License, or # 10 # (at your option) any later version. # 11 # # 12 # This program is distributed in the hope that it will be useful, # 13 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 15 # GNU General Public License for more details. # 16 # # 17 # You should have received a copy of the GNU General Public License # 18 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 19 # # 20 ############################################################################### 21 22 # Module docstring. 23 """Module for handling all types of structural superimpositions.""" 24 25 # Python module imports. 26 from copy import deepcopy 27 from math import pi 28 from numpy import diag, dot, eye, float64, outer, sign, transpose, zeros 29 from numpy.linalg import det, norm, svd 30 31 # relax module import. 32 from lib.structure.statistics import calc_mean_structure 33 from lib.geometry.rotations import R_to_axis_angle, R_to_euler_zyz 34 35

36 -def find_centroid(coords):

37 """Calculate the centroid of the structural coordinates. 38 39 @keyword coord: The atomic coordinates. 40 @type coord: numpy rank-2, Nx3 array 41 @return: The centroid. 42 @rtype: numpy rank-1, 3D array 43 """ 44 45 # The sum. 46 centroid = coords.sum(0) / coords.shape[0] 47 48 # Return. 49 return centroid

50 51

52 -def fit_to_first(models=None, coord=None, centroid=None):

53 """Superimpose a set of structural models using the fit to first algorithm. 54 55 @keyword models: The list of models to superimpose. 56 @type models: list of int 57 @keyword coord: The list of coordinates of all models to superimpose. The first index is the models, the second is the atomic positions, and the third is the xyz coordinates. 58 @type coord: list of numpy rank-2, Nx3 arrays 59 @keyword centroid: An alternative position of the centroid to allow for different superpositions, for example of pivot point motions. 60 @type centroid: list of float or numpy rank-1, 3D array 61 @return: The lists of translation vectors, rotation matrices, and rotation pivots. 62 @rtype: list of numpy rank-1 3D arrays, list of numpy rank-2 3D arrays, list of numpy rank-1 3D arrays 63 """ 64 65 # Print out. 66 print("\nSuperimposition of structural models %s using the 'fit to first' algorithm." % models) 67 68 # Init (there is no transformation for the first model). 69 T_list = [zeros(3, float64)] 70 R_list = [eye(3, dtype=float64)] 71 pivot_list = [zeros(3, float64)] 72 73 # Loop over the ending models. 74 for i in range(1, len(models)): 75 # Calculate the displacements (Kabsch algorithm). 76 trans_vect, trans_dist, R, axis, angle, pivot = kabsch(name_from='model %s'%models[0], name_to='model %s'%models[i], coord_from=coord[i], coord_to=coord[0], centroid=centroid) 77 78 # Store the transforms. 79 T_list.append(trans_vect) 80 R_list.append(R) 81 pivot_list.append(pivot) 82 83 # Return the transform data. 84 return T_list, R_list, pivot_list

85 86

87 -def fit_to_mean(models=None, coord=None, centroid=None, verbosity=1):

88 """Superimpose a set of structural models using the fit to first algorithm. 89 90 @keyword models: The list of models to superimpose. 91 @type models: list of int 92 @keyword coord: The list of coordinates of all models to superimpose. The first index is the models, the second is the atomic positions, and the third is the xyz coordinates. 93 @type coord: list of numpy rank-2, Nx3 arrays 94 @keyword centroid: An alternative position of the centroid to allow for different superpositions, for example of pivot point motions. 95 @type centroid: list of float or numpy rank-1, 3D array 96 @keyword verbosity: The amount of information to print out. If 0, nothing will be printed. 97 @type verbosity: int 98 @return: The lists of translation vectors, rotation matrices, and rotation pivots. 99 @rtype: list of numpy rank-1 3D arrays, list of numpy rank-2 3D arrays, list of numpy rank-1 3D arrays 100 """ 101 102 # Print out. 103 if verbosity: 104 print("\nSuperimposition of structural models %s using the 'fit to mean' algorithm." % models) 105 106 # Duplicate the coordinates. 107 orig_coord = deepcopy(coord) 108 109 # Initialise the displacement lists. 110 T_list = [] 111 R_list = [] 112 pivot_list = [] 113 114 # Initialise the mean structure. 115 N = len(coord[0]) 116 mean = zeros((N, 3), float64) 117 118 # Iterative fitting to mean. 119 converged = False 120 iter = 0 121 while not converged: 122 # Print out. 123 if verbosity: 124 print("\nIteration %i of the algorithm." % iter) 125 print("%-10s%-25s%-25s" % ("Model", "Translation (Angstrom)", "Rotation (deg)")) 126 127 # Calculate the mean structure. 128 calc_mean_structure(coord, mean) 129 130 # Fit each model to the mean. 131 converged = True 132 for i in range(len(models)): 133 # Calculate the displacements (Kabsch algorithm). 134 trans_vect, trans_dist, R, axis, angle, pivot = kabsch(name_from='model %s'%models[0], name_to='mean', coord_from=coord[i], coord_to=mean, centroid=centroid, verbosity=0) 135 136 # Table printout. 137 if verbosity: 138 print("%-10i%25.3g%25.3g" % (i, trans_dist, (angle / 2.0 / pi * 360.0))) 139 140 # Shift the coordinates. 141 for j in range(N): 142 # Translate. 143 coord[i][j] += trans_vect 144 145 # The pivot to atom vector. 146 coord[i][j] -= pivot 147 148 # Rotation. 149 coord[i][j] = dot(R, coord[i][j]) 150 151 # The new position. 152 coord[i][j] += pivot 153 154 # Convergence test. 155 if trans_dist > 1e-10 or angle > 1e-10: 156 converged = False 157 158 # Increment the iteration number. 159 iter += 1 160 161 # Perform the fit once from the original coordinates to obtain the full transforms. 162 for i in range(len(models)): 163 # Calculate the displacements (Kabsch algorithm). 164 trans_vect, trans_dist, R, axis, angle, pivot = kabsch(name_from='model %s'%models[i], name_to='the mean structure', coord_from=orig_coord[i], coord_to=mean, centroid=centroid, verbosity=0) 165 166 # Store the transforms. 167 T_list.append(trans_vect) 168 R_list.append(R) 169 pivot_list.append(pivot) 170 171 # Return the transform data. 172 return T_list, R_list, pivot_list

173 174

175 -def kabsch(name_from=None, name_to=None, coord_from=None, coord_to=None, centroid=None, verbosity=1):

176 """Calculate the rotational and translational displacements between the two given coordinate sets. 177 178 This uses the U{Kabsch algorithm<http://en.wikipedia.org/wiki/Kabsch_algorithm>}. 179 180 181 @keyword name_from: The name of the starting structure, used for the printouts. 182 @type name_from: str 183 @keyword name_to: The name of the ending structure, used for the printouts. 184 @type name_to: str 185 @keyword coord_from: The list of atomic coordinates for the starting structure. 186 @type coord_from: numpy rank-2, Nx3 array 187 @keyword coord_to: The list of atomic coordinates for the ending structure. 188 @type coord_to: numpy rank-2, Nx3 array 189 @keyword centroid: An alternative position of the centroid, used for studying pivoted systems. 190 @type centroid: list of float or numpy rank-1, 3D array 191 @return: The translation vector T, translation distance d, rotation matrix R, rotation axis r, rotation angle theta, and the rotational pivot defined as the centroid of the ending structure. 192 @rtype: numpy rank-1 3D array, float, numpy rank-2 3D array, numpy rank-1 3D array, float, numpy rank-1 3D array 193 """ 194 195 # Calculate the centroids. 196 if centroid != None: 197 centroid_from = centroid 198 centroid_to = centroid 199 else: 200 centroid_from = find_centroid(coord_from) 201 centroid_to = find_centroid(coord_to) 202 203 # The translation. 204 trans_vect = centroid_to - centroid_from 205 trans_dist = norm(trans_vect) 206 207 # Calculate the rotation. 208 R = kabsch_rotation(coord_from=coord_from, coord_to=coord_to, centroid_from=centroid_from, centroid_to=centroid_to) 209 axis, angle = R_to_axis_angle(R) 210 a, b, g = R_to_euler_zyz(R) 211 212 # Print out. 213 if verbosity >= 1: 214 print("\n\nCalculating the rotational and translational displacements from %s to %s using the Kabsch algorithm.\n" % (name_from, name_to)) 215 print("Start centroid: [%20.15f, %20.15f, %20.15f]" % (centroid_from[0], centroid_from[1], centroid_from[2])) 216 print("End centroid: [%20.15f, %20.15f, %20.15f]" % (centroid_to[0], centroid_to[1], centroid_to[2])) 217 print("Translation vector: [%20.15f, %20.15f, %20.15f]" % (trans_vect[0], trans_vect[1], trans_vect[2])) 218 print("Translation distance: %.15f" % trans_dist) 219 print("Rotation matrix:") 220 print(" [[%20.15f, %20.15f, %20.15f]" % (R[0, 0], R[0, 1], R[0, 2])) 221 print(" [%20.15f, %20.15f, %20.15f]" % (R[1, 0], R[1, 1], R[1, 2])) 222 print(" [%20.15f, %20.15f, %20.15f]]" % (R[2, 0], R[2, 1], R[2, 2])) 223 print("Rotation axis: [%20.15f, %20.15f, %20.15f]" % (axis[0], axis[1], axis[2])) 224 print("Rotation euler angles: [%20.15f, %20.15f, %20.15f]" % (a, b, g)) 225 print("Rotation angle (deg): %.15f" % (angle / 2.0 / pi * 360.0)) 226 227 # Return the data. 228 return trans_vect, trans_dist, R, axis, angle, centroid_to

229 230

231 -def kabsch_rotation(coord_from=None, coord_to=None, centroid_from=None, centroid_to=None):

232 """Calculate the rotation via SVD. 233 234 @keyword coord_from: The list of atomic coordinates for the starting structure. 235 @type coord_from: numpy rank-2, Nx3 array 236 @keyword coord_to: The list of atomic coordinates for the ending structure. 237 @type coord_to: numpy rank-2, Nx3 array 238 @keyword centroid_from: The starting centroid. 239 @type centroid_from: numpy rank-1, 3D array 240 @keyword centroid_to: The ending centroid. 241 @type centroid_to: numpy rank-1, 3D array 242 @return: The rotation matrix. 243 @rtype: numpy rank-2, 3D array 244 """ 245 246 # Initialise the covariance matrix A. 247 A = zeros((3, 3), float64) 248 249 # Loop over the atoms. 250 for i in range(coord_from.shape[0]): 251 # The positions shifted to the origin. 252 orig_from = coord_from[i] - centroid_from 253 orig_to = coord_to[i] - centroid_to 254 255 # The outer product. 256 A += outer(orig_from, orig_to) 257 258 # SVD. 259 U, S, V = svd(A) 260 261 # The handedness of the covariance matrix. 262 d = sign(det(A)) 263 D = diag([1, 1, d]) 264 265 # The rotation. 266 R = dot(transpose(V), dot(D, transpose(U))) 267 268 # Return the rotation. 269 return R

270

Source Code for Module lib.structure.superimpose