specific_analyses.frame

1 ############################################################################### 2 # # 3 # Copyright (C) 2009-2011,2013-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 all of the frame order specific user functions.""" 24 25 # Python module imports. 26 from copy import deepcopy 27 from math import pi 28 from numpy import argsort, array, float64, ones, transpose, zeros 29 from warnings import warn 30 31 # relax module imports. 32 from lib.arg_check import is_float_array 33 from lib.check_types import is_float 34 from lib.errors import RelaxError, RelaxFault 35 from lib.frame_order.simulation import brownian, mode_distribution, uniform_distribution 36 from lib.frame_order.variables import MODEL_DOUBLE_ROTOR, MODEL_ISO_CONE, MODEL_LIST, MODEL_LIST_FREE_ROTORS, MODEL_LIST_ISO_CONE, MODEL_LIST_PSEUDO_ELLIPSE, MODEL_LIST_RESTRICTED_TORSION, MODEL_PSEUDO_ELLIPSE, MODEL_RIGID 37 from lib.geometry.coord_transform import cartesian_to_spherical 38 from lib.geometry.rotations import R_to_euler_zyz 39 from lib.io import open_write_file 40 from lib.warnings import RelaxWarning 41 from pipe_control.pipes import check_pipe 42 from specific_analyses.frame_order.checks import check_domain, check_model, check_parameters, check_pivot 43 from specific_analyses.frame_order.data import domain_moving, generate_pivot 44 from specific_analyses.frame_order.geometric import average_position, create_ave_pos, create_geometric_rep, generate_axis_system 45 from specific_analyses.frame_order.optimisation import count_sobol_points 46 from specific_analyses.frame_order.parameters import assemble_param_vector, update_model 47 48

49 -def decompose(root="decomposed", dir=None, atom_id=None, model=1, total=100, reverse=False, mirror=False, force=True):

50 """Structural representation of the individual frame order motional components. 51 52 @keyword root: The file root for the PDB files created. Each motional component will be represented by a different PDB file appended with '_mode1.pdb', '_mode2.pdb', '_mode3.pdb', etc. 53 @type root: str 54 @keyword dir: The directory name to place the file into. 55 @type dir: str or None 56 @keyword atom_id: The atom identification string to allow the decomposition to be applied to subset of all atoms. 57 @type atom_id: None or str 58 @keyword model: Only one model from an analysed ensemble of structures can be used for the decomposition, as the corresponding PDB file consists of one model per state. 59 @type model: int 60 @keyword total: The total number of structures to distribute along the motional modes. This overrides a default fixed angle incrementation. 61 @type total: int 62 @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. 63 @type reverse: bool or list of bool 64 @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. 65 @type mirror: bool 66 @keyword force: A flag which, if set to True, will overwrite the any pre-existing file. 67 @type force: bool 68 """ 69 70 # Printout. 71 print("PDB representation of the individual components of the frame order motions.") 72 73 # Checks. 74 check_pipe() 75 check_model() 76 check_domain() 77 check_parameters() 78 check_pivot() 79 80 # Skip any unsupported models. 81 unsupported = [MODEL_RIGID, MODEL_DOUBLE_ROTOR] 82 if cdp.model in unsupported: 83 print("Skipping the unsupported '%s' model." % cdp.model) 84 return 85 86 # Initialise the angle vector (cone opening angle 1, cone opening angle 2, torsion angle). 87 angles = zeros(3, float64) 88 89 # Cone opening. 90 if cdp.model in MODEL_LIST_ISO_CONE: 91 angles[0] = angles[1] = cdp.cone_theta 92 elif cdp.model in MODEL_LIST_PSEUDO_ELLIPSE: 93 angles[0] = cdp.cone_theta_y 94 angles[1] = cdp.cone_theta_x 95 96 # Non-zero torsion angle. 97 if cdp.model in MODEL_LIST_FREE_ROTORS: 98 angles[2] = pi 99 elif cdp.model in MODEL_LIST_RESTRICTED_TORSION: 100 angles[2] = cdp.cone_sigma_max 101 102 # The motional eigenframe. 103 frame = generate_axis_system() 104 105 # Mode ordering from largest to smallest. 106 indices = argsort(angles) 107 angles = angles[indices[::-1]] 108 frame = transpose(transpose(frame)[indices[::-1]]) 109 110 # The pivot point. 111 pivot = generate_pivot(order=1, pdb_limit=True) 112 113 # Expand the list of reversal flags to all 3 modes, if required. 114 if isinstance(reverse, list): 115 reverse_list = reverse 116 for i in range(3 - len(reverse)): 117 reverse_list.append(False) 118 else: 119 reverse_list = [] 120 for i in range(3): 121 reverse_list.append(reverse) 122 123 # Loop over each mode. 124 for i in range(3): 125 # Skip modes with no motion. 126 if angles[i] < 1e-7: 127 continue 128 129 # Open the output file. 130 file_name = "%s_mode%i.pdb" % (root, i+1) 131 file = open_write_file(file_name=file_name, dir=dir, force=force) 132 133 # The structure. 134 structure = deepcopy(cdp.structure) 135 if structure.num_models() > 1: 136 structure.collapse_ensemble(model_num=model) 137 138 # Shift to the average position. 139 average_position(structure=structure, models=[None]) 140 141 # Create the representation. 142 mode_distribution(file=file, structure=structure, axis=frame[:,i], angle=angles[i], pivot=pivot, atom_id=domain_moving(), total=total, reverse=reverse_list[i], mirror=mirror) 143 144 # Close the file. 145 file.close()

146 147

148 -def distribute(file="distribution.pdb.bz2", dir=None, atom_id=None, total=1000, max_rotations=100000, model=1, force=True):

149 """Create a uniform distribution of structures for the frame order motions. 150 151 @keyword file: The PDB file for storing the frame order motional distribution. The compression is determined automatically by the file extensions '*.pdb', '*.pdb.gz', and '*.pdb.bz2'. 152 @type file: str 153 @keyword dir: The directory name to place the file into. 154 @type dir: str or None 155 @keyword atom_id: The atom identification string to allow the distribution to be a subset of all atoms. 156 @type atom_id: None or str 157 @keyword total: The total number of states/model/structures in the distribution. 158 @type total: int 159 @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. 160 @type max_rotations: int 161 @keyword model: Only one model from an analysed ensemble of structures can be used for the distribution, as the corresponding PDB file consists of one model per state. 162 @type model: int 163 @keyword force: A flag which, if set to True, will overwrite the any pre-existing file. 164 @type force: bool 165 """ 166 167 # Printout. 168 print("Uniform distribution of structures representing the frame order motions.") 169 170 # Check the total. 171 if total > 9999: 172 raise RelaxError("A maximum of 9999 models is allowed in the PDB format.") 173 174 # Checks. 175 check_pipe() 176 check_model() 177 check_domain() 178 check_parameters() 179 check_pivot() 180 181 # Skip the rigid model. 182 if cdp.model == MODEL_RIGID: 183 print("Skipping the rigid model.") 184 return 185 186 # Open the output file. 187 file = open_write_file(file_name=file, dir=dir, force=force) 188 189 # The parameter values. 190 values = assemble_param_vector() 191 params = {} 192 i = 0 193 for name in cdp.params: 194 params[name] = values[i] 195 i += 1 196 197 # The structure. 198 structure = deepcopy(cdp.structure) 199 if structure.num_models() > 1: 200 structure.collapse_ensemble(model_num=model) 201 202 # The pivot points. 203 num_states = 1 204 if cdp.model == MODEL_DOUBLE_ROTOR: 205 num_states = 2 206 pivot = zeros((num_states, 3), float64) 207 for i in range(num_states): 208 pivot[i] = generate_pivot(order=i+1, pdb_limit=True) 209 210 # Shift to the average position. 211 average_position(structure=structure, models=[None]) 212 213 # The motional eigenframe. 214 frame = generate_axis_system() 215 216 # Only work with a subset. 217 if atom_id: 218 # The inverted selection. 219 selection = structure.selection(atom_id=atom_id, inv=True) 220 221 # Delete the data. 222 structure.delete(selection=selection, verbosity=0) 223 224 # Create the distribution. 225 uniform_distribution(file=file, model=cdp.model, structure=structure, parameters=params, eigenframe=frame, pivot=pivot, atom_id=domain_moving(), total=total, max_rotations=max_rotations) 226 227 # Close the file. 228 file.close()

229 230

231 -def pdb_model(ave_pos="ave_pos", rep="frame_order", dir=None, compress_type=0, size=30.0, inc=36, model=1, force=False):

232 """Create 3 different PDB files for representing the frame order dynamics of the system. 233 234 @keyword ave_pos: The file root for the average molecule structure. 235 @type ave_pos: str or None 236 @keyword rep: The file root of the PDB representation of the frame order dynamics to create. 237 @type rep: str or None 238 @keyword dist: The file root which will contain multiple models spanning the full dynamics distribution of the frame order model. 239 @type dist: str or None 240 @keyword dir: The name of the directory to place the PDB file into. 241 @type dir: str 242 @keyword compress_type: The compression type. The integer values correspond to the compression type: 0, no compression; 1, Bzip2 compression; 2, Gzip compression. 243 @type compress_type: int 244 @keyword size: The size of the geometric object in Angstroms. 245 @type size: float 246 @keyword inc: The number of increments for the filling of the cone objects. 247 @type inc: int 248 @keyword model: Only one model from an analysed ensemble can be used for the PDB representation of the Monte Carlo simulations, as these consists of one model per simulation. 249 @type model: int 250 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten. 251 @type force: bool 252 """ 253 254 # Check that at least one PDB file name is given. 255 if not ave_pos and not rep and not dist: 256 raise RelaxError("Minimally one PDB file name must be supplied.") 257 258 # Test if the current data pipe exists. 259 check_pipe() 260 261 # Create the average position structure. 262 if ave_pos: 263 create_ave_pos(file=ave_pos, dir=dir, compress_type=compress_type, model=model, force=force) 264 265 # Nothing more to do for the rigid model. 266 if cdp.model == MODEL_RIGID: 267 return 268 269 # Create the geometric representation. 270 if rep: 271 create_geometric_rep(file=rep, dir=dir, compress_type=compress_type, size=size, inc=inc, force=force)

272 273

274 -def permute_axes(permutation='A'):

275 """Permute the axes of the motional eigenframe to switch between local minima. 276 277 @keyword permutation: The permutation to use. This can be either 'A' or 'B' to select between the 3 permutations, excluding the current combination. 278 @type permutation: str 279 """ 280 281 # Check that the model is valid. 282 allowed = MODEL_LIST_ISO_CONE + MODEL_LIST_PSEUDO_ELLIPSE 283 if cdp.model not in allowed: 284 raise RelaxError("The permutation of the motional eigenframe is only valid for the frame order models %s." % allowed) 285 286 # Check that the model parameters are setup. 287 if cdp.model in MODEL_LIST_ISO_CONE: 288 if not hasattr(cdp, 'cone_theta') or not is_float(cdp.cone_theta): 289 raise RelaxError("The parameter values are not set up.") 290 else: 291 if not hasattr(cdp, 'cone_theta_y') or not is_float(cdp.cone_theta_y): 292 raise RelaxError("The parameter values are not set up.") 293 294 # The iso cones only have one permutation. 295 if cdp.model in MODEL_LIST_ISO_CONE and permutation == 'B': 296 raise RelaxError("The isotropic cones only have one permutation.") 297 298 # The angles. 299 cone_sigma_max = 0.0 300 if cdp.model in MODEL_LIST_RESTRICTED_TORSION: 301 cone_sigma_max = cdp.cone_sigma_max 302 elif cdp.model in MODEL_LIST_FREE_ROTORS: 303 cone_sigma_max = pi 304 if cdp.model in MODEL_LIST_ISO_CONE: 305 angles = array([cdp.cone_theta, cdp.cone_theta, cone_sigma_max], float64) 306 else: 307 angles = array([cdp.cone_theta_x, cdp.cone_theta_y, cone_sigma_max], float64) 308 x, y, z = angles 309 310 # The axis system. 311 axes = generate_axis_system() 312 313 # Start printout for the isotropic cones. 314 if cdp.model in MODEL_LIST_ISO_CONE: 315 print("\nOriginal parameters:") 316 print("%-20s %20.10f" % ("cone_theta", cdp.cone_theta)) 317 print("%-20s %20.10f" % ("cone_sigma_max", cone_sigma_max)) 318 print("%-20s %20.10f" % ("axis_theta", cdp.axis_theta)) 319 print("%-20s %20.10f" % ("axis_phi", cdp.axis_phi)) 320 print("%-20s\n%s" % ("cone axis", axes[:, 2])) 321 print("%-20s\n%s" % ("full axis system", axes)) 322 print("\nPermutation '%s':" % permutation) 323 324 # Start printout for the pseudo-ellipses. 325 else: 326 print("\nOriginal parameters:") 327 print("%-20s %20.10f" % ("cone_theta_x", cdp.cone_theta_x)) 328 print("%-20s %20.10f" % ("cone_theta_y", cdp.cone_theta_y)) 329 print("%-20s %20.10f" % ("cone_sigma_max", cone_sigma_max)) 330 print("%-20s %20.10f" % ("eigen_alpha", cdp.eigen_alpha)) 331 print("%-20s %20.10f" % ("eigen_beta", cdp.eigen_beta)) 332 print("%-20s %20.10f" % ("eigen_gamma", cdp.eigen_gamma)) 333 print("%-20s\n%s" % ("eigenframe", axes)) 334 print("\nPermutation '%s':" % permutation) 335 336 # The axis inversion structure. 337 inv = ones(3, float64) 338 339 # The starting condition x <= y <= z. 340 if x <= y and y <= z: 341 # Printout. 342 print("%-20s %-20s" % ("Starting condition", "x <= y <= z")) 343 344 # The cone angle and axes permutations. 345 if permutation == 'A': 346 perm_angles = [0, 2, 1] 347 perm_axes = [2, 1, 0] 348 inv[perm_axes[2]] = -1.0 349 else: 350 perm_angles = [1, 2, 0] 351 perm_axes = [2, 0, 1] 352 353 # The starting condition x <= z <= y. 354 elif x <= z and z <= y: 355 # Printout. 356 print("%-20s %-20s" % ("Starting condition", "x <= z <= y")) 357 358 # The cone angle and axes permutations. 359 if permutation == 'A': 360 perm_angles = [0, 2, 1] 361 perm_axes = [2, 1, 0] 362 inv[perm_axes[2]] = -1.0 363 else: 364 perm_angles = [2, 1, 0] 365 perm_axes = [0, 2, 1] 366 inv[perm_axes[2]] = -1.0 367 368 # The starting condition z <= x <= y. 369 elif z <= x and x <= y: 370 # Printout. 371 print("%-20s %-20s" % ("Starting condition", "z <= x <= y")) 372 373 # The cone angle and axes permutations. 374 if permutation == 'A': 375 perm_angles = [2, 0, 1] 376 perm_axes = [1, 2, 0] 377 else: 378 perm_angles = [2, 1, 0] 379 perm_axes = [0, 2, 1] 380 inv[perm_axes[2]] = -1.0 381 382 # Cannot be here. 383 else: 384 raise RelaxFault 385 386 # Printout. 387 print("%-20s %-20s" % ("Cone angle permutation", perm_angles)) 388 print("%-20s %-20s" % ("Axes permutation", perm_axes)) 389 390 # Permute the angles. 391 if cdp.model in MODEL_LIST_ISO_CONE: 392 cdp.cone_theta = (angles[perm_angles[0]] + angles[perm_angles[1]]) / 2.0 393 else: 394 cdp.cone_theta_x = angles[perm_angles[0]] 395 cdp.cone_theta_y = angles[perm_angles[1]] 396 if cdp.model in MODEL_LIST_RESTRICTED_TORSION: 397 cdp.cone_sigma_max = angles[perm_angles[2]] 398 elif cdp.model in MODEL_LIST_FREE_ROTORS: 399 cdp.cone_sigma_max = pi 400 401 # Permute the axes (iso cone). 402 if cdp.model in MODEL_LIST_ISO_CONE: 403 # Convert the y-axis to spherical coordinates (the x-axis would be ok too, or any vector in the x-y plane due to symmetry of the original permutation). 404 axis_new = axes[:, 1] 405 r, cdp.axis_theta, cdp.axis_phi = cartesian_to_spherical(axis_new) 406 407 # Permute the axes (pseudo-ellipses). 408 else: 409 axes_new = transpose(array([inv[0]*axes[:, perm_axes[0]], inv[1]*axes[:, perm_axes[1]], inv[2]*axes[:, perm_axes[2]]], float64)) 410 411 # Convert the permuted frame to Euler angles and store them. 412 cdp.eigen_alpha, cdp.eigen_beta, cdp.eigen_gamma = R_to_euler_zyz(axes_new) 413 414 # End printout. 415 if cdp.model in MODEL_LIST_ISO_CONE: 416 print("\nPermuted parameters:") 417 print("%-20s %20.10f" % ("cone_theta", cdp.cone_theta)) 418 if cdp.model == MODEL_ISO_CONE: 419 print("%-20s %20.10f" % ("cone_sigma_max", cdp.cone_sigma_max)) 420 print("%-20s %20.10f" % ("axis_theta", cdp.axis_theta)) 421 print("%-20s %20.10f" % ("axis_phi", cdp.axis_phi)) 422 print("%-20s\n%s" % ("cone axis", axis_new)) 423 else: 424 print("\nPermuted parameters:") 425 print("%-20s %20.10f" % ("cone_theta_x", cdp.cone_theta_x)) 426 print("%-20s %20.10f" % ("cone_theta_y", cdp.cone_theta_y)) 427 if cdp.model == MODEL_PSEUDO_ELLIPSE: 428 print("%-20s %20.10f" % ("cone_sigma_max", cdp.cone_sigma_max)) 429 print("%-20s %20.10f" % ("eigen_alpha", cdp.eigen_alpha)) 430 print("%-20s %20.10f" % ("eigen_beta", cdp.eigen_beta)) 431 print("%-20s %20.10f" % ("eigen_gamma", cdp.eigen_gamma)) 432 print("%-20s\n%s" % ("eigenframe", axes_new))

433 434

435 -def pivot(pivot=None, order=1, fix=False):

436 """Set the pivot point for the 2 body motion. 437 438 @keyword pivot: The pivot point of the two bodies (domains, etc.) in the structural coordinate system. 439 @type pivot: list of num 440 @keyword order: The ordinal number of the pivot point. The value of 1 is for the first pivot point, the value of 2 for the second pivot point, and so on. 441 @type order: int 442 @keyword fix: A flag specifying if the pivot point should be fixed during optimisation. 443 @type fix: bool 444 """ 445 446 # Test if the current data pipe exists. 447 check_pipe() 448 449 # Store the fixed flag. 450 cdp.pivot_fixed = fix 451 452 # No pivot given, so update the model if needed and quit. 453 if pivot is None: 454 if hasattr(cdp, 'model'): 455 update_model() 456 return 457 458 # Convert the pivot to a numpy array. 459 pivot = array(pivot, float64) 460 461 # Check the pivot validity. 462 is_float_array(pivot, name='pivot point', size=3) 463 464 # Store the pivot point and fixed flag. 465 if order == 1: 466 cdp.pivot_x = pivot[0] 467 cdp.pivot_y = pivot[1] 468 cdp.pivot_z = pivot[2] 469 else: 470 # The variable names. 471 name_x = 'pivot_x_%i' % order 472 name_y = 'pivot_y_%i' % order 473 name_z = 'pivot_z_%i' % order 474 475 # Store the variables. 476 setattr(cdp, name_x, pivot[0]) 477 setattr(cdp, name_y, pivot[1]) 478 setattr(cdp, name_z, pivot[2]) 479 480 # Update the model. 481 if hasattr(cdp, 'model'): 482 update_model()

483 484

485 -def quad_int(flag=False):

486 """Turn the high precision Scipy quadratic numerical integration on or off. 487 488 @keyword flag: The flag which if True will perform high precision numerical integration via the scipy.integrate quad(), dblquad() and tplquad() integration methods rather than the rough quasi-random numerical integration. 489 @type flag: bool 490 """ 491 492 # Test if the current data pipe exists. 493 check_pipe() 494 495 # Store the flag. 496 cdp.quad_int = flag

497 498

499 -def ref_domain(ref=None):

500 """Set the reference domain for the frame order, multi-domain models. 501 502 @param ref: The reference domain. 503 @type ref: str 504 """ 505 506 # Checks. 507 check_pipe() 508 check_domain(domain=ref, escalate=0) 509 510 # Set the reference domain. 511 cdp.ref_domain = ref 512 513 # Update the model. 514 if hasattr(cdp, 'model'): 515 update_model()

516 517

518 -def select_model(model=None):

519 """Select the Frame Order model. 520 521 @param model: The Frame Order model. This can be one of 'pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor', 'iso cone', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor', 'rigid', 'free rotor', 'double rotor'. 522 @type model: str 523 """ 524 525 # Test if the current data pipe exists. 526 check_pipe() 527 528 # Test if the model name exists. 529 if not model in MODEL_LIST: 530 raise RelaxError("The model name '%s' is invalid, it must be one of %s." % (model, MODEL_LIST)) 531 532 # Set the model 533 cdp.model = model 534 535 # Initialise the list of model parameters. 536 cdp.params = [] 537 538 # Set the integration method if needed. 539 if not hasattr(cdp, 'quad_int'): 540 # Scipy quadratic numerical integration. 541 if cdp.model in []: 542 cdp.quad_int = True 543 544 # Quasi-random numerical integration. 545 else: 546 cdp.quad_int = False 547 548 # Update the model. 549 update_model()

550 551

552 -def simulate(file="simulation.pdb.bz2", dir=None, step_size=2.0, snapshot=10, total=1000, model=1, force=True):

553 """Pseudo-Brownian dynamics simulation of the frame order motions. 554 555 @keyword file: The PDB file for storing the frame order pseudo-Brownian dynamics simulation. The compression is determined automatically by the file extensions '*.pdb', '*.pdb.gz', and '*.pdb.bz2'. 556 @type file: str 557 @keyword dir: The directory name to place the file into. 558 @type dir: str or None 559 @keyword step_size: The rotation will be of a random direction but with this fixed angle. The value is in degrees. 560 @type step_size: float 561 @keyword snapshot: The number of steps in the simulation when snapshots will be taken. 562 @type snapshot: int 563 @keyword total: The total number of snapshots to take before stopping the simulation. 564 @type total: int 565 @keyword model: Only one model from an analysed ensemble of structures can be used for the pseudo-Brownian simulation, as the simulation and corresponding PDB file consists of one model per simulation. 566 @type model: int 567 @keyword force: A flag which, if set to True, will overwrite the any pre-existing file. 568 @type force: bool 569 """ 570 571 # Printout. 572 print("Pseudo-Brownian dynamics simulation of the frame order motions.") 573 574 # Checks. 575 check_pipe() 576 check_model() 577 check_domain() 578 check_parameters() 579 check_pivot() 580 581 # Skip the rigid model. 582 if cdp.model == MODEL_RIGID: 583 print("Skipping the rigid model.") 584 return 585 586 # Open the output file. 587 file = open_write_file(file_name=file, dir=dir, force=force) 588 589 # The parameter values. 590 values = assemble_param_vector() 591 params = {} 592 i = 0 593 for name in cdp.params: 594 params[name] = values[i] 595 i += 1 596 597 # The structure. 598 structure = deepcopy(cdp.structure) 599 if structure.num_models() > 1: 600 structure.collapse_ensemble(model_num=model) 601 602 # The pivot points. 603 num_states = 1 604 if cdp.model == MODEL_DOUBLE_ROTOR: 605 num_states = 2 606 pivot = zeros((num_states, 3), float64) 607 for i in range(num_states): 608 pivot[i] = generate_pivot(order=i+1, pdb_limit=True) 609 610 # Shift to the average position. 611 average_position(structure=structure, models=[None]) 612 613 # The motional eigenframe. 614 frame = generate_axis_system() 615 616 # Create the distribution. 617 brownian(file=file, model=cdp.model, structure=structure, parameters=params, eigenframe=frame, pivot=pivot, atom_id=domain_moving(), step_size=step_size, snapshot=snapshot, total=total) 618 619 # Close the file. 620 file.close()

621 622

623 -def sobol_setup(max_num=200, oversample=100):

624 """Oversampling setup for the quasi-random Sobol' sequence used for numerical PCS integration. 625 626 @keyword max_num: The maximum number of integration points N. 627 @type max_num: int 628 @keyword oversample: The oversampling factor Ov used for the N * Ov * 10**M, where M is the number of dimensions or torsion-tilt angles for the system. 629 @type oversample: int 630 """ 631 632 # Test if the current data pipe exists. 633 check_pipe() 634 635 # Throw a warning to the user if not enough points are being used. 636 if max_num < 200: 637 warn(RelaxWarning("To obtain reliable results in a frame order analysis, the maximum number of integration points should be greater than 200.")) 638 639 # Store the values. 640 cdp.sobol_max_points = max_num 641 cdp.sobol_oversample = oversample 642 643 # Count the number of Sobol' points for the current model. 644 count_sobol_points()

645

Source Code for Module specific_analyses.frame_order.uf