specific_analyses.frame

1 ############################################################################### 2 # # 3 # Copyright (C) 2009-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 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 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, 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 distribute(file="distribution.pdb.bz2", dir=None, atom_id=None, total=1000, max_rotations=100000, model=1, force=True):

50 """Create a uniform distribution of structures for the frame order motions. 51 52 @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'. 53 @type file: 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 distribution to be a subset of all atoms. 57 @type atom_id: None or str 58 @keyword total: The total number of states/model/structures in the distribution. 59 @type total: int 60 @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. 61 @type max_rotations: int 62 @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. 63 @type model: int 64 @keyword force: A flag which, if set to True, will overwrite the any pre-existing file. 65 @type force: bool 66 """ 67 68 # Printout. 69 print("Uniform distribution of structures representing the frame order motions.") 70 71 # Check the total. 72 if total > 9999: 73 raise RelaxError("A maximum of 9999 models is allowed in the PDB format.") 74 75 # Checks. 76 check_pipe() 77 check_model() 78 check_domain() 79 check_parameters() 80 check_pivot() 81 82 # Skip the rigid model. 83 if cdp.model == MODEL_RIGID: 84 print("Skipping the rigid model.") 85 return 86 87 # Open the output file. 88 file = open_write_file(file_name=file, dir=dir, force=force) 89 90 # The parameter values. 91 values = assemble_param_vector() 92 params = {} 93 i = 0 94 for name in cdp.params: 95 params[name] = values[i] 96 i += 1 97 98 # The structure. 99 structure = deepcopy(cdp.structure) 100 if structure.num_models() > 1: 101 structure.collapse_ensemble(model_num=model) 102 103 # The pivot points. 104 num_states = 1 105 if cdp.model == MODEL_DOUBLE_ROTOR: 106 num_states = 2 107 pivot = zeros((num_states, 3), float64) 108 for i in range(num_states): 109 pivot[i] = generate_pivot(order=i+1, pdb_limit=True) 110 111 # Shift to the average position. 112 average_position(structure=structure, models=[None]) 113 114 # The motional eigenframe. 115 frame = generate_axis_system() 116 117 # Only work with a subset. 118 if atom_id: 119 # The inverted selection. 120 selection = structure.selection(atom_id=atom_id, inv=True) 121 122 # Delete the data. 123 structure.delete(selection=selection, verbosity=0) 124 125 # Create the distribution. 126 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) 127 128 # Close the file. 129 file.close()

130 131

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

133 """Create 3 different PDB files for representing the frame order dynamics of the system. 134 135 @keyword ave_pos: The file root for the average molecule structure. 136 @type ave_pos: str or None 137 @keyword rep: The file root of the PDB representation of the frame order dynamics to create. 138 @type rep: str or None 139 @keyword dist: The file root which will contain multiple models spanning the full dynamics distribution of the frame order model. 140 @type dist: str or None 141 @keyword dir: The name of the directory to place the PDB file into. 142 @type dir: str 143 @keyword compress_type: The compression type. The integer values correspond to the compression type: 0, no compression; 1, Bzip2 compression; 2, Gzip compression. 144 @type compress_type: int 145 @keyword size: The size of the geometric object in Angstroms. 146 @type size: float 147 @keyword inc: The number of increments for the filling of the cone objects. 148 @type inc: int 149 @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. 150 @type model: int 151 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten. 152 @type force: bool 153 """ 154 155 # Check that at least one PDB file name is given. 156 if not ave_pos and not rep and not dist: 157 raise RelaxError("Minimally one PDB file name must be supplied.") 158 159 # Test if the current data pipe exists. 160 check_pipe() 161 162 # Create the average position structure. 163 if ave_pos: 164 create_ave_pos(file=ave_pos, dir=dir, compress_type=compress_type, model=model, force=force) 165 166 # Nothing more to do for the rigid model. 167 if cdp.model == MODEL_RIGID: 168 return 169 170 # Create the geometric representation. 171 if rep: 172 create_geometric_rep(file=rep, dir=dir, compress_type=compress_type, size=size, inc=inc, force=force)

173 174

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

176 """Permute the axes of the motional eigenframe to switch between local minima. 177 178 @keyword permutation: The permutation to use. This can be either 'A' or 'B' to select between the 3 permutations, excluding the current combination. 179 @type permutation: str 180 """ 181 182 # Check that the model is valid. 183 allowed = MODEL_LIST_ISO_CONE + MODEL_LIST_PSEUDO_ELLIPSE 184 if cdp.model not in allowed: 185 raise RelaxError("The permutation of the motional eigenframe is only valid for the frame order models %s." % allowed) 186 187 # Check that the model parameters are setup. 188 if cdp.model in MODEL_LIST_ISO_CONE: 189 if not hasattr(cdp, 'cone_theta') or not is_float(cdp.cone_theta): 190 raise RelaxError("The parameter values are not set up.") 191 else: 192 if not hasattr(cdp, 'cone_theta_y') or not is_float(cdp.cone_theta_y): 193 raise RelaxError("The parameter values are not set up.") 194 195 # The iso cones only have one permutation. 196 if cdp.model in MODEL_LIST_ISO_CONE and permutation == 'B': 197 raise RelaxError("The isotropic cones only have one permutation.") 198 199 # The angles. 200 cone_sigma_max = 0.0 201 if cdp.model in MODEL_LIST_RESTRICTED_TORSION: 202 cone_sigma_max = cdp.cone_sigma_max 203 elif cdp.model in MODEL_LIST_FREE_ROTORS: 204 cone_sigma_max = pi 205 if cdp.model in MODEL_LIST_ISO_CONE: 206 angles = array([cdp.cone_theta, cdp.cone_theta, cone_sigma_max], float64) 207 else: 208 angles = array([cdp.cone_theta_x, cdp.cone_theta_y, cone_sigma_max], float64) 209 x, y, z = angles 210 211 # The axis system. 212 axes = generate_axis_system() 213 214 # Start printout for the isotropic cones. 215 if cdp.model in MODEL_LIST_ISO_CONE: 216 print("\nOriginal parameters:") 217 print("%-20s %20.10f" % ("cone_theta", cdp.cone_theta)) 218 print("%-20s %20.10f" % ("cone_sigma_max", cone_sigma_max)) 219 print("%-20s %20.10f" % ("axis_theta", cdp.axis_theta)) 220 print("%-20s %20.10f" % ("axis_phi", cdp.axis_phi)) 221 print("%-20s\n%s" % ("cone axis", axes[:, 2])) 222 print("%-20s\n%s" % ("full axis system", axes)) 223 print("\nPermutation '%s':" % permutation) 224 225 # Start printout for the pseudo-ellipses. 226 else: 227 print("\nOriginal parameters:") 228 print("%-20s %20.10f" % ("cone_theta_x", cdp.cone_theta_x)) 229 print("%-20s %20.10f" % ("cone_theta_y", cdp.cone_theta_y)) 230 print("%-20s %20.10f" % ("cone_sigma_max", cone_sigma_max)) 231 print("%-20s %20.10f" % ("eigen_alpha", cdp.eigen_alpha)) 232 print("%-20s %20.10f" % ("eigen_beta", cdp.eigen_beta)) 233 print("%-20s %20.10f" % ("eigen_gamma", cdp.eigen_gamma)) 234 print("%-20s\n%s" % ("eigenframe", axes)) 235 print("\nPermutation '%s':" % permutation) 236 237 # The axis inversion structure. 238 inv = ones(3, float64) 239 240 # The starting condition x <= y <= z. 241 if x <= y and y <= z: 242 # Printout. 243 print("%-20s %-20s" % ("Starting condition", "x <= y <= z")) 244 245 # The cone angle and axes permutations. 246 if permutation == 'A': 247 perm_angles = [0, 2, 1] 248 perm_axes = [2, 1, 0] 249 inv[perm_axes[2]] = -1.0 250 else: 251 perm_angles = [1, 2, 0] 252 perm_axes = [2, 0, 1] 253 254 # The starting condition x <= z <= y. 255 elif x <= z and z <= y: 256 # Printout. 257 print("%-20s %-20s" % ("Starting condition", "x <= z <= y")) 258 259 # The cone angle and axes permutations. 260 if permutation == 'A': 261 perm_angles = [0, 2, 1] 262 perm_axes = [2, 1, 0] 263 inv[perm_axes[2]] = -1.0 264 else: 265 perm_angles = [2, 1, 0] 266 perm_axes = [0, 2, 1] 267 inv[perm_axes[2]] = -1.0 268 269 # The starting condition z <= x <= y. 270 elif z <= x and x <= y: 271 # Printout. 272 print("%-20s %-20s" % ("Starting condition", "z <= x <= y")) 273 274 # The cone angle and axes permutations. 275 if permutation == 'A': 276 perm_angles = [2, 0, 1] 277 perm_axes = [1, 2, 0] 278 else: 279 perm_angles = [2, 1, 0] 280 perm_axes = [0, 2, 1] 281 inv[perm_axes[2]] = -1.0 282 283 # Cannot be here. 284 else: 285 raise RelaxFault 286 287 # Printout. 288 print("%-20s %-20s" % ("Cone angle permutation", perm_angles)) 289 print("%-20s %-20s" % ("Axes permutation", perm_axes)) 290 291 # Permute the angles. 292 if cdp.model in MODEL_LIST_ISO_CONE: 293 cdp.cone_theta = (angles[perm_angles[0]] + angles[perm_angles[1]]) / 2.0 294 else: 295 cdp.cone_theta_x = angles[perm_angles[0]] 296 cdp.cone_theta_y = angles[perm_angles[1]] 297 if cdp.model in MODEL_LIST_RESTRICTED_TORSION: 298 cdp.cone_sigma_max = angles[perm_angles[2]] 299 elif cdp.model in MODEL_LIST_FREE_ROTORS: 300 cdp.cone_sigma_max = pi 301 302 # Permute the axes (iso cone). 303 if cdp.model in MODEL_LIST_ISO_CONE: 304 # 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). 305 axis_new = axes[:, 1] 306 r, cdp.axis_theta, cdp.axis_phi = cartesian_to_spherical(axis_new) 307 308 # Permute the axes (pseudo-ellipses). 309 else: 310 axes_new = transpose(array([inv[0]*axes[:, perm_axes[0]], inv[1]*axes[:, perm_axes[1]], inv[2]*axes[:, perm_axes[2]]], float64)) 311 312 # Convert the permuted frame to Euler angles and store them. 313 cdp.eigen_alpha, cdp.eigen_beta, cdp.eigen_gamma = R_to_euler_zyz(axes_new) 314 315 # End printout. 316 if cdp.model in MODEL_LIST_ISO_CONE: 317 print("\nPermuted parameters:") 318 print("%-20s %20.10f" % ("cone_theta", cdp.cone_theta)) 319 if cdp.model == MODEL_ISO_CONE: 320 print("%-20s %20.10f" % ("cone_sigma_max", cdp.cone_sigma_max)) 321 print("%-20s %20.10f" % ("axis_theta", cdp.axis_theta)) 322 print("%-20s %20.10f" % ("axis_phi", cdp.axis_phi)) 323 print("%-20s\n%s" % ("cone axis", axis_new)) 324 else: 325 print("\nPermuted parameters:") 326 print("%-20s %20.10f" % ("cone_theta_x", cdp.cone_theta_x)) 327 print("%-20s %20.10f" % ("cone_theta_y", cdp.cone_theta_y)) 328 if cdp.model == MODEL_PSEUDO_ELLIPSE: 329 print("%-20s %20.10f" % ("cone_sigma_max", cdp.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_new))

334 335

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

337 """Set the pivot point for the 2 body motion. 338 339 @keyword pivot: The pivot point of the two bodies (domains, etc.) in the structural coordinate system. 340 @type pivot: list of num 341 @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. 342 @type order: int 343 @keyword fix: A flag specifying if the pivot point should be fixed during optimisation. 344 @type fix: bool 345 """ 346 347 # Test if the current data pipe exists. 348 check_pipe() 349 350 # Store the fixed flag. 351 cdp.pivot_fixed = fix 352 353 # No pivot given, so update the model if needed and quit. 354 if pivot is None: 355 if hasattr(cdp, 'model'): 356 update_model() 357 return 358 359 # Convert the pivot to a numpy array. 360 pivot = array(pivot, float64) 361 362 # Check the pivot validity. 363 is_float_array(pivot, name='pivot point', size=3) 364 365 # Store the pivot point and fixed flag. 366 if order == 1: 367 cdp.pivot_x = pivot[0] 368 cdp.pivot_y = pivot[1] 369 cdp.pivot_z = pivot[2] 370 else: 371 # The variable names. 372 name_x = 'pivot_x_%i' % order 373 name_y = 'pivot_y_%i' % order 374 name_z = 'pivot_z_%i' % order 375 376 # Store the variables. 377 setattr(cdp, name_x, pivot[0]) 378 setattr(cdp, name_y, pivot[1]) 379 setattr(cdp, name_z, pivot[2]) 380 381 # Update the model. 382 if hasattr(cdp, 'model'): 383 update_model()

384 385

386 -def quad_int(flag=False):

387 """Turn the high precision Scipy quadratic numerical integration on or off. 388 389 @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. 390 @type flag: bool 391 """ 392 393 # Test if the current data pipe exists. 394 check_pipe() 395 396 # Store the flag. 397 cdp.quad_int = flag

398 399

400 -def ref_domain(ref=None):

401 """Set the reference domain for the frame order, multi-domain models. 402 403 @param ref: The reference domain. 404 @type ref: str 405 """ 406 407 # Checks. 408 check_pipe() 409 check_domain(domain=ref, escalate=0) 410 411 # Set the reference domain. 412 cdp.ref_domain = ref 413 414 # Update the model. 415 if hasattr(cdp, 'model'): 416 update_model()

417 418

419 -def select_model(model=None):

420 """Select the Frame Order model. 421 422 @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'. 423 @type model: str 424 """ 425 426 # Test if the current data pipe exists. 427 check_pipe() 428 429 # Test if the model name exists. 430 if not model in MODEL_LIST: 431 raise RelaxError("The model name '%s' is invalid, it must be one of %s." % (model, MODEL_LIST)) 432 433 # Set the model 434 cdp.model = model 435 436 # Initialise the list of model parameters. 437 cdp.params = [] 438 439 # Set the integration method if needed. 440 if not hasattr(cdp, 'quad_int'): 441 # Scipy quadratic numerical integration. 442 if cdp.model in []: 443 cdp.quad_int = True 444 445 # Quasi-random numerical integration. 446 else: 447 cdp.quad_int = False 448 449 # Update the model. 450 update_model()

451 452

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

454 """Pseudo-Brownian dynamics simulation of the frame order motions. 455 456 @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'. 457 @type file: str 458 @keyword dir: The directory name to place the file into. 459 @type dir: str or None 460 @keyword step_size: The rotation will be of a random direction but with this fixed angle. The value is in degrees. 461 @type step_size: float 462 @keyword snapshot: The number of steps in the simulation when snapshots will be taken. 463 @type snapshot: int 464 @keyword total: The total number of snapshots to take before stopping the simulation. 465 @type total: int 466 @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. 467 @type model: int 468 @keyword force: A flag which, if set to True, will overwrite the any pre-existing file. 469 @type force: bool 470 """ 471 472 # Printout. 473 print("Pseudo-Brownian dynamics simulation of the frame order motions.") 474 475 # Checks. 476 check_pipe() 477 check_model() 478 check_domain() 479 check_parameters() 480 check_pivot() 481 482 # Skip the rigid model. 483 if cdp.model == MODEL_RIGID: 484 print("Skipping the rigid model.") 485 return 486 487 # Open the output file. 488 file = open_write_file(file_name=file, dir=dir, force=force) 489 490 # The parameter values. 491 values = assemble_param_vector() 492 params = {} 493 i = 0 494 for name in cdp.params: 495 params[name] = values[i] 496 i += 1 497 498 # The structure. 499 structure = deepcopy(cdp.structure) 500 if structure.num_models() > 1: 501 structure.collapse_ensemble(model_num=model) 502 503 # The pivot points. 504 num_states = 1 505 if cdp.model == MODEL_DOUBLE_ROTOR: 506 num_states = 2 507 pivot = zeros((num_states, 3), float64) 508 for i in range(num_states): 509 pivot[i] = generate_pivot(order=i+1, pdb_limit=True) 510 511 # Shift to the average position. 512 average_position(structure=structure, models=[None]) 513 514 # The motional eigenframe. 515 frame = generate_axis_system() 516 517 # Create the distribution. 518 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) 519 520 # Close the file. 521 file.close()

522 523

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

525 """Oversampling setup for the quasi-random Sobol' sequence used for numerical PCS integration. 526 527 @keyword max_num: The maximum number of integration points N. 528 @type max_num: int 529 @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. 530 @type oversample: int 531 """ 532 533 # Test if the current data pipe exists. 534 check_pipe() 535 536 # Throw a warning to the user if not enough points are being used. 537 if max_num < 200: 538 warn(RelaxWarning("To obtain reliable results in a frame order analysis, the maximum number of integration points should be greater than 200.")) 539 540 # Store the values. 541 cdp.sobol_max_points = max_num 542 cdp.sobol_oversample = oversample 543 544 # Count the number of Sobol' points for the current model. 545 count_sobol_points()

546

Source Code for Module specific_analyses.frame_order.uf