generic_fns.structure.superimpose

1 ############################################################################### 2 # # 3 # Copyright (C) 2011 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, sqrt, transpose, zeros 29 from numpy.linalg import det, norm, svd 30 31 # relax module import. 32 from maths_fns.rotation_matrix import R_to_axis_angle 33 34

35 -class Pivot_finder:

36 """Class for finding the optimal pivot point for motions between the given models.""" 37

38 - def __init__(self, models, coord):

39 """Set up the class for pivot point optimisation. 40 41 @keyword models: The list of models to use. If set to None, then all models will be used. 42 @type models: list of int or None 43 @keyword coord: The array of molecular coordinates. The first dimension corresponds to the model, the second the atom, the third the coordinate. 44 @type coord: rank-3 numpy array 45 """ 46 47 # Store the args. 48 self.models = models 49 self.coord = coord 50 51 # Store a copy of the coordinates for restoration. 52 self.orig_coord = deepcopy(coord)

53 54

55 - def func(self, params):

56 """Target function for the optimisation of the motional pivot point. 57 58 @param params: The parameter vector from the optimisation algorithm. 59 @type params: list 60 @return: The target function value defined as the combined RMSD value. 61 @rtype: float 62 """ 63 64 # The fit to mean algorithm. 65 T, R, pivot = fit_to_mean(models=self.models, coord=self.coord, centroid=params, verbosity=0) 66 67 # The RMSD. 68 val = atomic_rmsd(self.coord) 69 70 # Restore the coordinates. 71 self.coord = deepcopy(self.orig_coord) 72 73 # Return the RMSD. 74 return val

75 76 77

78 -def atomic_rmsd(coord):

79 """Determine the RMSD for the given atomic coordinates. 80 81 This is the per atom RMSD to the mean structure. 82 83 84 @keyword coord: The array of molecular coordinates. The first dimension corresponds to the model, the second the atom, the third the coordinate. 85 @type coord: rank-3 numpy array 86 """ 87 88 # Init. 89 M = len(coord) 90 N = len(coord[0]) 91 model_rmsd = zeros(M, float64) 92 mean = zeros((N, 3), float64) 93 94 # Calculate the mean structure. 95 calc_mean_structure(coord, mean) 96 97 # Loop over the models. 98 for i in range(M): 99 # Loop over the atoms. 100 for j in range(N): 101 # The vector connecting the mean to model atom. 102 vect = mean[j] - coord[i][j] 103 104 # The atomic RMSD. 105 model_rmsd[i] += norm(vect)**2 106 107 # Normalise, and sqrt. 108 model_rmsd[i] = sqrt(model_rmsd[i] / N) 109 110 # Return the average RMSD. 111 return sum(model_rmsd) / M

112 113

114 -def calc_mean_structure(coord=None, mean=None):

115 """Average the coordinates. 116 117 @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. 118 @type coord: list of numpy rank-2, Nx3 arrays 119 @keyword mean: The data storage for the mean structure. 120 @type mean: numpy rank-2, Nx3 array 121 """ 122 123 # The number of atoms. 124 N = len(coord[0]) 125 M = len(coord) 126 127 # Clear the mean data structure. 128 for i in range(N): 129 mean[i] = [0.0, 0.0, 0.0] 130 131 # Loop over the atoms. 132 for i in range(N): 133 # Loop over the models. 134 for j in range(M): 135 mean[i] += coord[j][i] 136 137 # Average. 138 mean[i] = mean[i] / M

139 140

141 -def find_centroid(coords):

142 """Calculate the centroid of the structural coordinates. 143 144 @keyword coord: The atomic coordinates. 145 @type coord: numpy rank-2, Nx3 array 146 @return: The centroid. 147 @rtype: numpy rank-1, 3D array 148 """ 149 150 # The sum. 151 centroid = coords.sum(0) / coords.shape[0] 152 153 # Return. 154 return centroid

155 156

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

158 """Superimpose a set of structural models using the fit to first algorithm. 159 160 @keyword models: The list of models to superimpose. 161 @type models: list of int 162 @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. 163 @type coord: list of numpy rank-2, Nx3 arrays 164 @keyword centroid: An alternative position of the centroid to allow for different superpositions, for example of pivot point motions. 165 @type centroid: list of float or numpy rank-1, 3D array 166 @return: The lists of translation vectors, rotation matrices, and rotation pivots. 167 @rtype: list of numpy rank-1 3D arrays, list of numpy rank-2 3D arrays, list of numpy rank-1 3D arrays 168 """ 169 170 # Print out. 171 print("\nSuperimposition of structural models %s using the 'fit to first' algorithm." % models) 172 173 # Init (there is no transformation for the first model). 174 T_list = [zeros(3, float64)] 175 R_list = [eye(3, dtype=float64)] 176 pivot_list = [zeros(3, float64)] 177 178 # Loop over the ending models. 179 for i in range(1, len(models)): 180 # Calculate the displacements (Kabsch algorithm). 181 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) 182 183 # Store the transforms. 184 T_list.append(trans_vect) 185 R_list.append(R) 186 pivot_list.append(pivot) 187 188 # Return the transform data. 189 return T_list, R_list, pivot_list

190 191

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

193 """Superimpose a set of structural models using the fit to first algorithm. 194 195 @keyword models: The list of models to superimpose. 196 @type models: list of int 197 @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. 198 @type coord: list of numpy rank-2, Nx3 arrays 199 @keyword centroid: An alternative position of the centroid to allow for different superpositions, for example of pivot point motions. 200 @type centroid: list of float or numpy rank-1, 3D array 201 @keyword verbosity: The amount of information to print out. If 0, nothing will be printed. 202 @type verbosity: int 203 @return: The lists of translation vectors, rotation matrices, and rotation pivots. 204 @rtype: list of numpy rank-1 3D arrays, list of numpy rank-2 3D arrays, list of numpy rank-1 3D arrays 205 """ 206 207 # Print out. 208 if verbosity: 209 print("\nSuperimposition of structural models %s using the 'fit to mean' algorithm." % models) 210 211 # Duplicate the coordinates. 212 orig_coord = deepcopy(coord) 213 214 # Initialise the displacement lists. 215 T_list = [] 216 R_list = [] 217 pivot_list = [] 218 219 # Initialise the mean structure. 220 N = len(coord[0]) 221 mean = zeros((N, 3), float64) 222 223 # Iterative fitting to mean. 224 converged = False 225 iter = 0 226 while not converged: 227 # Print out. 228 if verbosity: 229 print("\nIteration %i of the algorithm." % iter) 230 print("%-10s%-25s%-25s" % ("Model", "Translation (Angstrom)", "Rotation (deg)")) 231 232 # Calculate the mean structure. 233 calc_mean_structure(coord, mean) 234 235 # Fit each model to the mean. 236 converged = True 237 for i in range(len(models)): 238 # Calculate the displacements (Kabsch algorithm). 239 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) 240 241 # Table printout. 242 if verbosity: 243 print("%-10i%25.3g%25.3g" % (i, trans_dist, (angle / 2.0 / pi * 360.0))) 244 245 # Shift the coordinates. 246 for j in range(N): 247 # Translate. 248 coord[i][j] += trans_vect 249 250 # The pivot to atom vector. 251 coord[i][j] -= pivot 252 253 # Rotation. 254 coord[i][j] = dot(R, coord[i][j]) 255 256 # The new position. 257 coord[i][j] += pivot 258 259 # Convergence test. 260 if trans_dist > 1e-10 or angle > 1e-10: 261 converged = False 262 263 # Increment the iteration number. 264 iter += 1 265 266 # Perform the fit once from the original coordinates to obtain the full transforms. 267 for i in range(len(models)): 268 # Calculate the displacements (Kabsch algorithm). 269 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) 270 271 # Store the transforms. 272 T_list.append(trans_vect) 273 R_list.append(R) 274 pivot_list.append(pivot) 275 276 # Return the transform data. 277 return T_list, R_list, pivot_list

278 279

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

281 """Calculate the rotational and translational displacements between the two given coordinate sets. 282 283 This uses the Kabsch algorithm (http://en.wikipedia.org/wiki/Kabsch_algorithm). 284 285 286 @keyword name_from: The name of the starting structure, used for the printouts. 287 @type name_from: str 288 @keyword name_to: The name of the ending structure, used for the printouts. 289 @type name_to: str 290 @keyword coord_from: The list of atomic coordinates for the starting structure. 291 @type coord_from: numpy rank-2, Nx3 array 292 @keyword coord_to: The list of atomic coordinates for the ending structure. 293 @type coord_to: numpy rank-2, Nx3 array 294 @keyword centroid: An alternative position of the centroid, used for studying pivoted systems. 295 @type centroid: list of float or numpy rank-1, 3D array 296 @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. 297 @rtype: numpy rank-1 3D array, float, numpy rank-2 3D array, numpy rank-1 3D array, float, numpy rank-1 3D array 298 """ 299 300 # Calculate the centroids. 301 if centroid != None: 302 centroid_from = centroid 303 centroid_to = centroid 304 else: 305 centroid_from = find_centroid(coord_from) 306 centroid_to = find_centroid(coord_to) 307 308 # The translation. 309 trans_vect = centroid_to - centroid_from 310 trans_dist = norm(trans_vect) 311 312 # Calculate the rotation. 313 R = kabsch_rotation(coord_from=coord_from, coord_to=coord_to, centroid_from=centroid_from, centroid_to=centroid_to) 314 axis, angle = R_to_axis_angle(R) 315 316 # Print out. 317 if verbosity >= 1: 318 print("\n\nCalculating the rotational and translational displacements from %s to %s using the Kabsch algorithm.\n" % (name_from, name_to)) 319 print("Start centroid: [%20.15f, %20.15f, %20.15f]" % (centroid_from[0], centroid_from[1], centroid_from[2])) 320 print("End centroid: [%20.15f, %20.15f, %20.15f]" % (centroid_to[0], centroid_to[1], centroid_to[2])) 321 print("Translation vector: [%20.15f, %20.15f, %20.15f]" % (trans_vect[0], trans_vect[1], trans_vect[2])) 322 print("Translation distance: %.15f" % trans_dist) 323 print("Rotation matrix:") 324 print(" [[%20.15f, %20.15f, %20.15f]" % (R[0, 0], R[0, 1], R[0, 2])) 325 print(" [%20.15f, %20.15f, %20.15f]" % (R[1, 0], R[1, 1], R[1, 2])) 326 print(" [%20.15f, %20.15f, %20.15f]]" % (R[2, 0], R[2, 1], R[2, 2])) 327 print("Rotation axis: [%20.15f, %20.15f, %20.15f]" % (axis[0], axis[1], axis[2])) 328 print("Rotation angle (deg): %.15f" % (angle / 2.0 / pi * 360.0)) 329 330 # Return the data. 331 return trans_vect, trans_dist, R, axis, angle, centroid_to

332 333

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

335 """Calculate the rotation via SVD. 336 337 @keyword coord_from: The list of atomic coordinates for the starting structure. 338 @type coord_from: numpy rank-2, Nx3 array 339 @keyword coord_to: The list of atomic coordinates for the ending structure. 340 @type coord_to: numpy rank-2, Nx3 array 341 @keyword centroid_from: The starting centroid. 342 @type centroid_from: numpy rank-1, 3D array 343 @keyword centroid_to: The ending centroid. 344 @type centroid_to: numpy rank-1, 3D array 345 @return: The rotation matrix. 346 @rtype: numpy rank-2, 3D array 347 """ 348 349 # Initialise the covariance matrix A. 350 A = zeros((3, 3), float64) 351 352 # Loop over the atoms. 353 for i in range(coord_from.shape[0]): 354 # The positions shifted to the origin. 355 orig_from = coord_from[i] - centroid_from 356 orig_to = coord_to[i] - centroid_to 357 358 # The outer product. 359 A += outer(orig_from, orig_to) 360 361 # SVD. 362 U, S, V = svd(A) 363 364 # The handedness of the covariance matrix. 365 d = sign(det(A)) 366 D = diag([1, 1, d]) 367 368 # The rotation. 369 R = dot(transpose(V), dot(D, transpose(U))) 370 371 # Return the rotation. 372 return R

373

Source Code for Module generic_fns.structure.superimpose