Package lib :: Package frame_order :: Module simulation
[hide private]
[frames] | no frames]

Source Code for Module lib.frame_order.simulation

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2014-2015,2018-2019 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 simulating the frame order motions.""" 
 24   
 25  # Python module imports. 
 26  from math import cos, modf, pi, sin, sqrt 
 27  from numpy import arange, array, concatenate, dot, eye, float64, linspace, transpose, zeros 
 28  import sys 
 29  from warnings import warn 
 30   
 31  # relax module imports. 
 32  from lib.errors import RelaxError 
 33  from lib.warnings import RelaxWarning 
 34  from lib.frame_order.variables import MODEL_DOUBLE_ROTOR 
 35  from lib.geometry.angles import wrap_angles 
 36  from lib.geometry.rotations import axis_angle_to_R, R_random_hypersphere, R_to_tilt_torsion, tilt_torsion_to_R 
 37  from lib.geometry.vectors import random_unit_vector 
 38   
 39   
40 -def brownian(file=None, model=None, structure=None, parameters={}, eigenframe=None, pivot=None, atom_id=None, step_size=2.0, snapshot=10, total=1000):
41 """Pseudo-Brownian dynamics simulation of the frame order motions. 42 43 @keyword file: The opened and writable file object to place the snapshots into. 44 @type file: str 45 @keyword structure: The internal structural object containing the domain to simulate as a single model. 46 @type structure: lib.structure.internal.object.Internal instance 47 @keyword model: The frame order model to simulate. 48 @type model: str 49 @keyword parameters: The dictionary of model parameter values. The key is the parameter name and the value is the value. 50 @type parameters: dict of float 51 @keyword eigenframe: The full 3D eigenframe of the frame order motions. 52 @type eigenframe: numpy rank-2, 3D float64 array 53 @keyword pivot: The list of pivot points of the frame order motions. 54 @type pivot: numpy rank-2 (N, 3) float64 array 55 @keyword atom_id: The atom ID string for the atoms in the structure to rotate - i.e. the moving domain. 56 @type atom_id: None or str 57 @keyword step_size: The rotation will be of a random direction but with this fixed angle. The value is in degrees. 58 @type step_size: float 59 @keyword snapshot: The number of steps in the simulation when snapshots will be taken. 60 @type snapshot: int 61 @keyword total: The total number of snapshots to take before stopping the simulation. 62 @type total: int 63 """ 64 65 # Check the structural object. 66 if structure.num_models() > 1: 67 raise RelaxError("Only a single model is supported.") 68 69 # Set the model number. 70 structure.set_model(model_orig=None, model_new=1) 71 72 # Generate the internal structural selection object. 73 selection = structure.selection(atom_id) 74 75 # The initial states and motional limits. 76 num_states = len(pivot) 77 states = zeros((num_states, 3, 3), float64) 78 theta_max = [] 79 sigma_max = [] 80 for i in range(num_states): 81 states[i] = eye(3) 82 theta_max.append(None) 83 sigma_max.append(None) 84 85 # Initialise the rotation matrix data structures. 86 vector = zeros(3, float64) 87 R = eye(3, dtype=float64) 88 step_size = step_size / 360.0 * 2.0 * pi 89 90 # Axis permutations. 91 perm = [None] 92 if model == MODEL_DOUBLE_ROTOR: 93 perm = [[2, 0, 1], [1, 2, 0]] 94 perm_rev = [[1, 2, 0], [2, 0, 1]] 95 96 # The maximum cone opening angles (isotropic cones). 97 if 'cone_theta' in parameters: 98 theta_max[0] = parameters['cone_theta'] 99 100 # The maximum cone opening angles (isotropic cones). 101 theta_x = None 102 theta_y = None 103 if 'cone_theta_x' in parameters: 104 theta_x = parameters['cone_theta_x'] 105 theta_y = parameters['cone_theta_y'] 106 107 # The maximum torsion angle. 108 if 'cone_sigma_max' in parameters: 109 sigma_max[0] = parameters['cone_sigma_max'] 110 elif 'free rotor' in model: 111 sigma_max[0] = pi 112 113 # The second torsion angle. 114 if 'cone_sigma_max_2' in parameters: 115 sigma_max[1] = parameters['cone_sigma_max_2'] 116 117 # Printout. 118 print("\nRunning the simulation:") 119 120 # Simulate. 121 current_snapshot = 1 122 step = 1 123 while True: 124 # End the simulation. 125 if current_snapshot == total: 126 print("\nEnd of simulation.") 127 break 128 129 # Loop over each state, or motional mode. 130 for i in range(num_states): 131 # The random vector. 132 random_unit_vector(vector) 133 134 # The rotation matrix. 135 axis_angle_to_R(vector, step_size, R) 136 137 # Shift the current state. 138 states[i] = dot(R, states[i]) 139 140 # Rotation in the eigenframe. 141 R_eigen = dot(transpose(eigenframe), dot(states[i], eigenframe)) 142 143 # Axis permutation to shift each rotation axis to Z. 144 if perm[i] != None: 145 for j in range(3): 146 R_eigen[:, j] = R_eigen[perm[i], j] 147 for j in range(3): 148 R_eigen[j, :] = R_eigen[j, perm[i]] 149 150 # The angles. 151 phi, theta, sigma = R_to_tilt_torsion(R_eigen) 152 sigma = wrap_angles(sigma, -pi, pi) 153 154 # Determine theta_max for the pseudo-ellipse models. 155 if theta_x != None: 156 theta_max[i] = 1.0 / sqrt((cos(phi) / theta_x)**2 + (sin(phi) / theta_y)**2) 157 158 # Set the cone opening angle to the maximum if outside of the limit. 159 if theta_max[i] != None: 160 if theta > theta_max[i]: 161 theta = theta_max[i] 162 163 # No tilt component. 164 else: 165 theta = 0.0 166 phi = 0.0 167 168 # Set the torsion angle to the maximum if outside of the limits. 169 if sigma_max[i] != None: 170 if sigma > sigma_max[i]: 171 sigma = sigma_max[i] 172 elif sigma < -sigma_max[i]: 173 sigma = -sigma_max[i] 174 else: 175 sigma = 0.0 176 177 # Reconstruct the rotation matrix, in the eigenframe, without sigma. 178 tilt_torsion_to_R(phi, theta, sigma, R_eigen) 179 180 # Reverse axis permutation to shift each rotation z-axis back. 181 if perm[i] != None: 182 for j in range(3): 183 R_eigen[:, j] = R_eigen[perm_rev[i], j] 184 for j in range(3): 185 R_eigen[j, :] = R_eigen[j, perm_rev[i]] 186 187 # Rotate back out of the eigenframe. 188 states[i] = dot(eigenframe, dot(R_eigen, transpose(eigenframe))) 189 190 # Take a snapshot. 191 if step == snapshot: 192 # Progress. 193 sys.stdout.write('.') 194 sys.stdout.flush() 195 196 # Increment the snapshot number. 197 current_snapshot += 1 198 199 # Copy the original structural data. 200 structure.add_model(model=current_snapshot, coords_from=1) 201 202 # Rotate the model. 203 for i in range(num_states): 204 structure.rotate(R=states[i], origin=pivot[i], model=current_snapshot, selection=selection) 205 206 # Reset the step counter. 207 step = 0 208 209 # Increment. 210 step += 1 211 212 # Save the result. 213 structure.write_pdb(file=file)
214 215
216 -def mode_distribution(file=None, structure=None, axis=None, angle=None, pivot=None, atom_id=None, angle_inc=2*pi/360, total=None, reverse=False, mirror=False):
217 """Linear distribution of a single component of the frame order motions. 218 219 @keyword file: The opened and wrlib/frame_order/simulation.pyitable file object to place the PDB models of the representation into. 220 @type file: str 221 @keyword structure: The internal structural object to convert into an ensemble along the mode of motion. 222 @type structure: lib.structure.internal.object.Internal instance 223 @keyword axis: The rotation axis. 224 @type axis: numpy 3D float64 array 225 @keyword angle: The rotation angle in radian (structures will be rotated +/- this angle). 226 @type angle: float 227 @keyword pivot: The pivot point for the given motional mode. 228 @type pivot: numpy 3D float64 array 229 @keyword atom_id: The atom ID string for the atoms in the structure to rotate - i.e. the moving domain. 230 @type atom_id: None or str 231 @keyword angle_inc: The angle between rotated representations. The default is 1 degree. 232 @type angle_inc: float 233 @keyword total: The total number of structures to distribute along the motional modes. This overrides angle_inc. 234 @type total: int 235 @keyword reverse: Set this to reverse the ordering of the models distributed along the motional mode. Use a list of Booleans to selectively reverse each motional mode. 236 @type reverse: bool or list of bool 237 @keyword mirror: Set this to have the models distributed along the motional mode shift from the negative angle to positive angle, and then return to the negative angle. 238 @type mirror: bool 239 """ 240 241 # Check the structural object. 242 if structure.num_models() > 1: 243 raise RelaxError("Only a single model is supported.") 244 245 # Set the model number. 246 structure.set_model(model_orig=None, model_new=1) 247 248 # Generate the internal structural selection object. 249 selection = structure.selection(atom_id) 250 251 # Initialise the rotation matrix data structures. 252 R = eye(3, dtype=float64) 253 254 # No angle handling. 255 if angle == 0.0 or total == 1: 256 angles = array([0.0]) 257 258 # Range of angles for a fixed number of structures. 259 elif total != None: 260 if mirror: 261 total = int(total / 2) + 1 262 angles = linspace(-angle, angle, total) 263 264 # Range of angles for a fixed angle between structures. 265 else: 266 angles = linspace(-angle, angle, int(angle/angle_inc)) 267 if not len(angles): 268 angles = array([-angle, 0.0, angle]) 269 270 # Angle reversal. 271 if reverse: 272 angles = angles[::-1] 273 274 # Angle mirroring. 275 if mirror: 276 angles2 = angles[:-1] 277 angles2 = angles2[::-1] 278 angles = concatenate((angles, angles2)) 279 280 # Generate the structures. 281 current_model = 1 282 for i in range(len(angles)): 283 # The rotation matrix. 284 axis_angle_to_R(axis, angles[i], R) 285 286 # Increment the snapshot number. 287 current_model += 1 288 289 # Copy the original structural data. 290 structure.add_model(model=current_model, coords_from=1) 291 292 # Rotate the model. 293 structure.rotate(R=R, origin=pivot, model=current_model, selection=selection) 294 295 # Delete the first model. 296 structure.delete(model=1) 297 298 # Save the result. 299 structure.write_pdb(file=file) 300 print("")
301 302
303 -def uniform_distribution(file=None, model=None, structure=None, parameters={}, eigenframe=None, pivot=None, atom_id=None, total=1000, max_rotations=100000):
304 """Uniform distribution of the frame order motions. 305 306 @keyword file: The opened and writable file object to place the PDB models of the distribution into. 307 @type file: str 308 @keyword structure: The internal structural object containing the domain to distribute as a single model. 309 @type structure: lib.structure.internal.object.Internal instance 310 @keyword model: The frame order model to distribute. 311 @type model: str 312 @keyword parameters: The dictionary of model parameter values. The key is the parameter name and the value is the value. 313 @type parameters: dict of float 314 @keyword eigenframe: The full 3D eigenframe of the frame order motions. 315 @type eigenframe: numpy rank-2, 3D float64 array 316 @keyword pivot: The list of pivot points of the frame order motions. 317 @type pivot: numpy rank-2 (N, 3) float64 array 318 @keyword atom_id: The atom ID string for the atoms in the structure to rotate - i.e. the moving domain. 319 @type atom_id: None or str 320 @keyword total: The total number of states in the distribution. 321 @type total: int 322 @keyword max_rotations: The maximum number of rotations to generate the distribution from. This prevents an execution for an infinite amount of time when a frame order amplitude parameter is close to zero so that the subset of all rotations within the distribution is close to zero. 323 @type max_rotations: int 324 """ 325 326 # Check the structural object. 327 if structure.num_models() > 1: 328 raise RelaxError("Only a single model is supported.") 329 330 # Set the model number. 331 structure.set_model(model_orig=None, model_new=1) 332 333 # Generate the internal structural selection object. 334 selection = structure.selection(atom_id) 335 336 # The initial states and motional limits. 337 num_states = len(pivot) 338 states = zeros((num_states, 3, 3), float64) 339 theta_max = [] 340 sigma_max = [] 341 for i in range(num_states): 342 states[i] = eye(3) 343 theta_max.append(None) 344 sigma_max.append(None) 345 346 # Initialise the rotation matrix data structures. 347 R = eye(3, dtype=float64) 348 349 # Axis permutations. 350 perm = [None] 351 if model == MODEL_DOUBLE_ROTOR: 352 perm = [[2, 0, 1], [1, 2, 0]] 353 perm_rev = [[1, 2, 0], [2, 0, 1]] 354 355 # The maximum cone opening angles (isotropic cones). 356 if 'cone_theta' in parameters: 357 theta_max[0] = parameters['cone_theta'] 358 359 # The maximum cone opening angles (isotropic cones). 360 theta_x = None 361 theta_y = None 362 if 'cone_theta_x' in parameters: 363 theta_x = parameters['cone_theta_x'] 364 theta_y = parameters['cone_theta_y'] 365 366 # The maximum torsion angle. 367 if 'cone_sigma_max' in parameters: 368 sigma_max[0] = parameters['cone_sigma_max'] 369 elif 'free rotor' in model: 370 sigma_max[0] = pi 371 372 # The second torsion angle. 373 if 'cone_sigma_max_2' in parameters: 374 sigma_max[1] = parameters['cone_sigma_max_2'] 375 376 # Printout. 377 print("\nGenerating the distribution:") 378 379 # Distribution. 380 current_state = 1 381 num = -1 382 while True: 383 # The total number of rotations. 384 num += 1 385 386 # End. 387 if current_state == total: 388 break 389 if num >= max_rotations: 390 sys.stdout.write('\n') 391 warn(RelaxWarning("Maximum number of rotations encountered - the distribution only contains %i states." % current_state)) 392 break 393 394 # Loop over each state, or motional mode. 395 inside = True 396 for i in range(num_states): 397 # The random rotation matrix. 398 R_random_hypersphere(R) 399 400 # Shift the current state. 401 states[i] = dot(R, states[i]) 402 403 # Rotation in the eigenframe. 404 R_eigen = dot(transpose(eigenframe), dot(states[i], eigenframe)) 405 406 # Axis permutation to shift each rotation axis to Z. 407 if perm[i] != None: 408 for j in range(3): 409 R_eigen[:, j] = R_eigen[perm[i], j] 410 for j in range(3): 411 R_eigen[j, :] = R_eigen[j, perm[i]] 412 413 # The angles. 414 phi, theta, sigma = R_to_tilt_torsion(R_eigen) 415 sigma = wrap_angles(sigma, -pi, pi) 416 417 # Determine theta_max for the pseudo-ellipse models. 418 if theta_x != None: 419 theta_max[i] = 1.0 / sqrt((cos(phi) / theta_x)**2 + (sin(phi) / theta_y)**2) 420 421 # The cone opening angle is outside of the limit. 422 if theta_max[i] != None: 423 if theta > theta_max[i]: 424 inside = False 425 426 # No tilt component. 427 else: 428 theta = 0.0 429 phi = 0.0 430 431 # The torsion angle is outside of the limits. 432 if sigma_max[i] != None: 433 if sigma > sigma_max[i]: 434 inside = False 435 elif sigma < -sigma_max[i]: 436 inside = False 437 else: 438 sigma = 0.0 439 440 # Reconstruct the rotation matrix, in the eigenframe, without sigma. 441 tilt_torsion_to_R(phi, theta, sigma, R_eigen) 442 443 # Reverse axis permutation to shift each rotation z-axis back. 444 if perm[i] != None: 445 for j in range(3): 446 R_eigen[:, j] = R_eigen[perm_rev[i], j] 447 for j in range(3): 448 R_eigen[j, :] = R_eigen[j, perm_rev[i]] 449 450 # Rotate back out of the eigenframe. 451 states[i] = dot(eigenframe, dot(R_eigen, transpose(eigenframe))) 452 453 # The state is outside of the distribution. 454 if not inside: 455 continue 456 457 # Progress. 458 sys.stdout.write('.') 459 sys.stdout.flush() 460 461 # Increment the snapshot number. 462 current_state += 1 463 464 # Copy the original structural data. 465 structure.add_model(model=current_state, coords_from=1) 466 467 # Rotate the model. 468 for i in range(num_states): 469 structure.rotate(R=states[i], origin=pivot[i], model=current_state, selection=selection) 470 471 # Save the result. 472 structure.write_pdb(file=file)
473