generic_fns.structure.superimpose

1 ############################################################################### 2 # # 3 # Copyright (C) 2011 Edward d'Auvergne # 4 # # 5 # This file is part of the program relax. # 6 # # 7 # relax 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 2 of the License, or # 10 # (at your option) any later version. # 11 # # 12 # relax 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 relax; if not, write to the Free Software # 19 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA # 20 # # 21 ############################################################################### 22 23 # Module docstring. 24 """Module for handling all types of structural superimpositions.""" 25 26 # Python module imports. 27 from copy import deepcopy 28 from math import pi 29 from numpy import diag, dot, eye, float64, outer, sign, sqrt, transpose, zeros 30 from numpy.linalg import det, norm, svd 31 32 # relax module import. 33 from maths_fns.rotation_matrix import R_to_axis_angle 34 35

36 -class Pivot_finder:

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

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

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

54 55

56 - def func(self, params):

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

76 77 78

79 -def atomic_rmsd(coord):

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

113 114

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

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

140 141

142 -def find_centroid(coords):

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

156 157

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

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

191 192

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

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

279 280

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

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

333 334

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

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

374

Source Code for Module generic_fns.structure.superimpose