lib.frame_order.simulation

1 ############################################################################### 2 # # 3 # Copyright (C) 2014-2015 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, pi, sin, sqrt 27 from numpy import dot, eye, float64, 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 uniform_distribution(file=None, model=None, structure=None, parameters={}, eigenframe=None, pivot=None, atom_id=None, total=1000, max_rotations=100000):

217 """Uniform distribution of the frame order motions. 218 219 @keyword file: The opened and writable file object to place the PDB models of the distribution into. 220 @type file: str 221 @keyword structure: The internal structural object containing the domain to distribute as a single model. 222 @type structure: lib.structure.internal.object.Internal instance 223 @keyword model: The frame order model to distribute. 224 @type model: str 225 @keyword parameters: The dictionary of model parameter values. The key is the parameter name and the value is the value. 226 @type parameters: dict of float 227 @keyword eigenframe: The full 3D eigenframe of the frame order motions. 228 @type eigenframe: numpy rank-2, 3D float64 array 229 @keyword pivot: The list of pivot points of the frame order motions. 230 @type pivot: numpy rank-2 (N, 3) float64 array 231 @keyword atom_id: The atom ID string for the atoms in the structure to rotate - i.e. the moving domain. 232 @type atom_id: None or str 233 @keyword total: The total number of states in the distribution. 234 @type total: int 235 @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. 236 @type max_rotations: int 237 """ 238 239 # Check the structural object. 240 if structure.num_models() > 1: 241 raise RelaxError("Only a single model is supported.") 242 243 # Set the model number. 244 structure.set_model(model_orig=None, model_new=1) 245 246 # Generate the internal structural selection object. 247 selection = structure.selection(atom_id) 248 249 # The initial states and motional limits. 250 num_states = len(pivot) 251 states = zeros((num_states, 3, 3), float64) 252 theta_max = [] 253 sigma_max = [] 254 for i in range(num_states): 255 states[i] = eye(3) 256 theta_max.append(None) 257 sigma_max.append(None) 258 259 # Initialise the rotation matrix data structures. 260 R = eye(3, dtype=float64) 261 262 # Axis permutations. 263 perm = [None] 264 if model == MODEL_DOUBLE_ROTOR: 265 perm = [[2, 0, 1], [1, 2, 0]] 266 perm_rev = [[1, 2, 0], [2, 0, 1]] 267 268 # The maximum cone opening angles (isotropic cones). 269 if 'cone_theta' in parameters: 270 theta_max[0] = parameters['cone_theta'] 271 272 # The maximum cone opening angles (isotropic cones). 273 theta_x = None 274 theta_y = None 275 if 'cone_theta_x' in parameters: 276 theta_x = parameters['cone_theta_x'] 277 theta_y = parameters['cone_theta_y'] 278 279 # The maximum torsion angle. 280 if 'cone_sigma_max' in parameters: 281 sigma_max[0] = parameters['cone_sigma_max'] 282 elif 'free rotor' in model: 283 sigma_max[0] = pi 284 285 # The second torsion angle. 286 if 'cone_sigma_max_2' in parameters: 287 sigma_max[1] = parameters['cone_sigma_max_2'] 288 289 # Printout. 290 print("\nGenerating the distribution:") 291 292 # Distribution. 293 current_state = 1 294 num = -1 295 while True: 296 # The total number of rotations. 297 num += 1 298 299 # End. 300 if current_state == total: 301 break 302 if num >= max_rotations: 303 sys.stdout.write('\n') 304 warn(RelaxWarning("Maximum number of rotations encountered - the distribution only contains %i states." % current_state)) 305 break 306 307 # Loop over each state, or motional mode. 308 inside = True 309 for i in range(num_states): 310 # The random rotation matrix. 311 R_random_hypersphere(R) 312 313 # Shift the current state. 314 states[i] = dot(R, states[i]) 315 316 # Rotation in the eigenframe. 317 R_eigen = dot(transpose(eigenframe), dot(states[i], eigenframe)) 318 319 # Axis permutation to shift each rotation axis to Z. 320 if perm[i] != None: 321 for j in range(3): 322 R_eigen[:, j] = R_eigen[perm[i], j] 323 for j in range(3): 324 R_eigen[j, :] = R_eigen[j, perm[i]] 325 326 # The angles. 327 phi, theta, sigma = R_to_tilt_torsion(R_eigen) 328 sigma = wrap_angles(sigma, -pi, pi) 329 330 # Determine theta_max for the pseudo-ellipse models. 331 if theta_x != None: 332 theta_max[i] = 1.0 / sqrt((cos(phi) / theta_x)**2 + (sin(phi) / theta_y)**2) 333 334 # The cone opening angle is outside of the limit. 335 if theta_max[i] != None: 336 if theta > theta_max[i]: 337 inside = False 338 339 # No tilt component. 340 else: 341 theta = 0.0 342 phi = 0.0 343 344 # The torsion angle is outside of the limits. 345 if sigma_max[i] != None: 346 if sigma > sigma_max[i]: 347 inside = False 348 elif sigma < -sigma_max[i]: 349 inside = False 350 else: 351 sigma = 0.0 352 353 # Reconstruct the rotation matrix, in the eigenframe, without sigma. 354 tilt_torsion_to_R(phi, theta, sigma, R_eigen) 355 356 # Reverse axis permutation to shift each rotation z-axis back. 357 if perm[i] != None: 358 for j in range(3): 359 R_eigen[:, j] = R_eigen[perm_rev[i], j] 360 for j in range(3): 361 R_eigen[j, :] = R_eigen[j, perm_rev[i]] 362 363 # Rotate back out of the eigenframe. 364 states[i] = dot(eigenframe, dot(R_eigen, transpose(eigenframe))) 365 366 # The state is outside of the distribution. 367 if not inside: 368 continue 369 370 # Progress. 371 sys.stdout.write('.') 372 sys.stdout.flush() 373 374 # Increment the snapshot number. 375 current_state += 1 376 377 # Copy the original structural data. 378 structure.add_model(model=current_state, coords_from=1) 379 380 # Rotate the model. 381 for i in range(num_states): 382 structure.rotate(R=states[i], origin=pivot[i], model=current_state, selection=selection) 383 384 # Save the result. 385 structure.write_pdb(file=file)

386

Source Code for Module lib.frame_order.simulation