specific_analyses.frame_order.geometric

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 handling the geometric representation of the frame order motions.""" 24 25 # Python module imports. 26 from copy import deepcopy 27 from math import pi 28 from numpy import cross, dot, eye, float64, zeros 29 from numpy.linalg import norm 30 import sys 31 from warnings import warn 32 33 # relax module imports. 34 from lib.errors import RelaxFault 35 from lib.frame_order.conversions import create_rotor_axis_alpha, create_rotor_axis_spherical 36 from lib.frame_order.variables import MODEL_DOUBLE_ROTOR, MODEL_FREE_ROTOR, MODEL_ISO_CONE, MODEL_ISO_CONE_FREE_ROTOR, MODEL_ISO_CONE_TORSIONLESS, MODEL_LIST_DOUBLE, MODEL_LIST_FREE_ROTORS, MODEL_LIST_ISO_CONE, MODEL_LIST_PSEUDO_ELLIPSE, MODEL_PSEUDO_ELLIPSE, MODEL_PSEUDO_ELLIPSE_FREE_ROTOR, MODEL_PSEUDO_ELLIPSE_TORSIONLESS, MODEL_ROTOR 37 from lib.geometry.rotations import euler_to_R_zyz 38 from lib.io import open_write_file 39 from lib.structure.cones import Iso_cone, Pseudo_elliptic 40 from lib.structure.geometric import generate_vector_residues 41 from lib.structure.internal.object import Internal 42 from lib.structure.represent.cone import cone 43 from lib.structure.represent.rotor import rotor 44 from lib.text.sectioning import subsection, subsubsection 45 from lib.warnings import RelaxWarning 46 from pipe_control.structure.mass import pipe_centre_of_mass 47 from specific_analyses.frame_order.checks import check_parameters 48 from specific_analyses.frame_order.data import domain_moving, generate_pivot 49 50

51 -def add_axes(structure=None, representation=None, size=None, sims=False):

52 """Add the axis system for the current frame order model to the structural object. 53 54 @keyword structure: The internal structural object to add the rotor objects to. 55 @type structure: lib.structure.internal.object.Internal instance 56 @keyword representation: The representation to create. If this is set to None or 'A', the standard representation will be created. If set to 'B', the axis system will be inverted. 57 @type representation: None or str 58 @keyword size: The size of the geometric object in Angstroms. 59 @type size: float 60 @keyword sims: A flag which if True will add the Monte Carlo simulation rotors to the structural object. There must be one model for each Monte Carlo simulation already present in the structural object. 61 @type sims: bool 62 """ 63 64 # Create the molecule. 65 mol_name = 'axes' 66 structure.add_molecule(name=mol_name) 67 68 # The transformation matrix (identity matrix or inversion matrix). 69 if representation == 'B': 70 T = -eye(3) 71 else: 72 T = eye(3) 73 74 # The models to loop over. 75 model_nums = [None] 76 sim_indices = [None] 77 if sims: 78 model_nums = [i+1 for i in range(cdp.sim_number)] 79 sim_indices = list(range(cdp.sim_number)) 80 81 # Loop over the models. 82 for i in range(len(model_nums)): 83 # Alias the molecule. 84 mol = structure.get_molecule(mol_name, model=model_nums[i]) 85 86 # The pivot points. 87 pivot1 = generate_pivot(order=1, sim_index=sim_indices[i], pdb_limit=True) 88 pivot2 = generate_pivot(order=2, sim_index=sim_indices[i], pdb_limit=True) 89 90 # A single z-axis, when no rotor object is present. 91 if cdp.model in [MODEL_ISO_CONE_TORSIONLESS]: 92 # Print out. 93 print("\nGenerating the z-axis system.") 94 95 # The axis. 96 if sims: 97 axis = create_rotor_axis_spherical(theta=cdp.axis_theta_sim[sim_indices[i]], phi=cdp.axis_phi_sim[sim_indices[i]]) 98 else: 99 axis = create_rotor_axis_spherical(theta=cdp.axis_theta, phi=cdp.axis_phi) 100 print(("Central axis: %s." % axis)) 101 102 # Transform the central axis. 103 axis = dot(T, axis) 104 105 # Generate the axis vectors. 106 print("\nGenerating the axis vectors.") 107 res_num = generate_vector_residues(mol=mol, vector=axis, atom_name='z-ax', res_name_vect='AXE', res_num=2, origin=pivot1, scale=size) 108 109 # The z-axis connecting two motional modes. 110 elif cdp.model in [MODEL_DOUBLE_ROTOR]: 111 # Printout. 112 print("\nGenerating the z-axis linking the two pivot points.") 113 114 # The axis. 115 axis = pivot1 - pivot2 116 print(("Interconnecting axis: %s." % axis)) 117 118 # Generate the axis vectors. 119 print("\nGenerating the axis vectors.") 120 res_num = generate_vector_residues(mol=mol, vector=axis, atom_name='z-ax', res_name_vect='AXE', res_num=1, origin=pivot2) 121 122 # The full axis system. 123 elif cdp.model in MODEL_LIST_PSEUDO_ELLIPSE: 124 # Print out. 125 print("\nGenerating the full axis system.") 126 127 # Add the pivot point. 128 mol.atom_add(atom_num=1, atom_name='R', res_name='AXE', res_num=1, pos=pivot1, element='C', pdb_record='HETATM') 129 130 # The axis system. 131 axes = generate_axis_system(sim_index=sim_indices[i]) 132 133 # Rotations and inversions. 134 axes = dot(T, axes) 135 136 # The axes to create. 137 label = ['x', 'y'] 138 if cdp.model in [MODEL_PSEUDO_ELLIPSE_TORSIONLESS]: 139 label = ['x', 'y', 'z'] 140 141 # Generate the axis vectors. 142 print("\nGenerating the axis vectors.") 143 for j in range(len(label)): 144 res_num = generate_vector_residues(mol=mol, vector=axes[:, j], atom_name='%s-ax'%label[j], res_name_vect='AXE', res_num=2, origin=pivot1, scale=size)

145 146

147 -def add_cones(structure=None, representation=None, size=None, inc=None, sims=False):

148 """Add the cone geometric object for the current frame order model to the structural object. 149 150 @keyword structure: The internal structural object to add the rotor objects to. 151 @type structure: lib.structure.internal.object.Internal instance 152 @keyword representation: The representation to create. If this is set to None or 'A', the standard representation will be created. If set to 'B', the axis system will be inverted. 153 @type representation: str 154 @keyword size: The size of the geometric object in Angstroms. 155 @type size: float 156 @keyword inc: The number of increments for the filling of the cone objects. 157 @type inc: int 158 @keyword sims: A flag which if True will add the Monte Carlo simulation pivots to the structural object. There must be one model for each Monte Carlo simulation already present in the structural object. 159 @type sims: bool 160 """ 161 162 # Create the molecule. 163 structure.add_molecule(name='cones') 164 165 # The transformation matrix (identity matrix or inversion matrix). 166 if representation == 'B': 167 T = -eye(3) 168 else: 169 T = eye(3) 170 171 # The models to loop over. 172 model_nums = [None] 173 sim_indices = [None] 174 if sims: 175 model_nums = [i+1 for i in range(cdp.sim_number)] 176 sim_indices = list(range(cdp.sim_number)) 177 178 # Loop over the models. 179 for i in range(len(model_nums)): 180 # Alias the molecule. 181 mol = structure.get_molecule('cones', model=model_nums[i]) 182 183 # The 1st pivot point. 184 pivot = generate_pivot(order=1, sim_index=sim_indices[i], pdb_limit=True) 185 186 # The rotation matrix (rotation from the z-axis to the cone axis). 187 R = generate_axis_system(sim_index=sim_indices[i]) 188 print(("Rotation matrix:\n%s" % R)) 189 190 # The transformation. 191 R = dot(T, R) 192 193 # The pseudo-ellipse cone object. 194 if cdp.model in MODEL_LIST_PSEUDO_ELLIPSE: 195 if sims: 196 cone_obj = Pseudo_elliptic(cdp.cone_theta_x_sim[sim_indices[i]], cdp.cone_theta_y_sim[sim_indices[i]]) 197 else: 198 cone_obj = Pseudo_elliptic(cdp.cone_theta_x, cdp.cone_theta_y) 199 200 # The isotropic cone object. 201 else: 202 # The angle. 203 if hasattr(cdp, 'cone_theta'): 204 if sims: 205 cone_theta = cdp.cone_theta_sim[sim_indices[i]] 206 else: 207 cone_theta = cdp.cone_theta 208 209 # The object. 210 cone_obj = Iso_cone(cone_theta) 211 212 # Create the cone. 213 cone(mol=mol, cone_obj=cone_obj, start_res=1, apex=pivot, R=R, inc=inc, scale=size, distribution='regular', axis_flag=False)

214 215

216 -def add_pivots(structure=None, sims=False):

217 """Add the pivots for the current frame order model to the structural object. 218 219 @keyword structure: The internal structural object to add the rotor objects to. 220 @type structure: lib.structure.internal.object.Internal instance 221 @keyword sims: A flag which if True will add the Monte Carlo simulation pivots to the structural object. There must be one model for each Monte Carlo simulation already present in the structural object. 222 @type sims: bool 223 """ 224 225 # Initialise. 226 pivots = [] 227 mols = [] 228 atom_names = [] 229 atom_nums = [] 230 231 # Create the molecule. 232 mol_name = 'pivots' 233 structure.add_molecule(name=mol_name) 234 235 # The models to loop over. 236 model_nums = [None] 237 sim_indices = [None] 238 if sims: 239 model_nums = [i+1 for i in range(cdp.sim_number)] 240 sim_indices = list(range(cdp.sim_number)) 241 242 # Loop over the models. 243 for i in range(len(model_nums)): 244 # Alias the molecule. 245 mol = structure.get_molecule(mol_name, model=model_nums[i]) 246 247 # The pivot points. 248 pivot1 = generate_pivot(order=1, sim_index=sim_indices[i], pdb_limit=True) 249 pivot2 = generate_pivot(order=2, sim_index=sim_indices[i], pdb_limit=True) 250 251 # Add the pivots for the double motion models. 252 if cdp.model in [MODEL_DOUBLE_ROTOR]: 253 # The 1st pivot. 254 mols.append(mol) 255 pivots.append(pivot1) 256 atom_names.append('Piv1') 257 atom_nums.append(1) 258 259 # The 2nd pivot. 260 mols.append(mol) 261 pivots.append(pivot2) 262 atom_names.append('Piv2') 263 atom_nums.append(2) 264 265 # Add the pivot for the single motion models. 266 else: 267 mols.append(mol) 268 pivots.append(pivot1) 269 atom_names.append('Piv') 270 atom_nums.append(1) 271 272 # Loop over the data, adding all pivots. 273 for i in range(len(mols)): 274 mols[i].atom_add(atom_num=atom_nums[i], atom_name=atom_names[i], res_name='PIV', res_num=1, pos=pivots[i], element='C', pdb_record='HETATM')

275 276

277 -def add_titles(structure=None, representation=None, displacement=40.0, sims=False):

278 """Add atoms to be used for the titles for the frame order geometric objects. 279 280 @keyword structure: The internal structural object to add the rotor objects to. 281 @type structure: lib.structure.internal.object.Internal instance 282 @keyword representation: The representation to title. 283 @type representation: None or str 284 @keyword displacement: The distance away from the pivot point, in Angstrom, to place the title. The simulation title will be shifted by a few extra Angstrom to avoid clashes. 285 @type displacement: float 286 @keyword sims: A flag which if True will add the Monte Carlo simulation pivots to the structural object. There must be one model for each Monte Carlo simulation already present in the structural object. 287 @type sims: bool 288 """ 289 290 # The title atom names. 291 atom_name = None 292 if representation == None and sims: 293 atom_name = 'mc' 294 elif representation == 'A': 295 if sims: 296 atom_name = 'mc-a' 297 else: 298 atom_name = 'a' 299 elif representation == 'B': 300 if sims: 301 atom_name = 'mc-b' 302 else: 303 atom_name = 'b' 304 305 # Nothing to do. 306 if atom_name == None: 307 return 308 309 # Avoid name overlaps. 310 if sims: 311 displacement += 3 312 313 # Create the molecule. 314 mol_name = 'titles' 315 structure.add_molecule(name=mol_name) 316 317 # The transformation matrix (identity matrix or inversion matrix). 318 if representation == 'B': 319 T = -eye(3) 320 else: 321 T = eye(3) 322 323 # The pivot points. 324 pivot1 = generate_pivot(order=1, pdb_limit=True) 325 pivot2 = generate_pivot(order=2, pdb_limit=True) 326 327 # The models to loop over. 328 model_nums = [None] 329 sim_indices = [None] 330 if sims: 331 model_nums = [i+1 for i in range(cdp.sim_number)] 332 sim_indices = list(range(cdp.sim_number)) 333 334 # Loop over the models. 335 for i in range(len(model_nums)): 336 # Alias the molecule. 337 mol = structure.get_molecule(mol_name, model=model_nums[i]) 338 339 # The axis system. 340 axes = generate_axis_system(sim_index=sim_indices[i]) 341 342 # Transform the central axis. 343 axis = dot(T, axes[:, 2]) 344 345 # The label position. 346 pos = pivot1 + displacement*axis 347 348 # Add the atom. 349 mol.atom_add(atom_name=atom_name, res_name='TLE', res_num=1, pos=pos, element='Ti', pdb_record='HETATM')

350 351

352 -def add_rotors(structure=None, representation=None, size=None, sims=False):

353 """Add all rotor objects for the current frame order model to the structural object. 354 355 @keyword structure: The internal structural object to add the rotor objects to. 356 @type structure: lib.structure.internal.object.Internal instance 357 @keyword representation: The representation to create. If this is set to None, the rotor will be dumbbell shaped with propellers at both ends. If set to 'A' or 'B', only half of the rotor will be shown. 358 @type representation: None or str 359 @keyword size: The size of the geometric object in Angstroms. 360 @type size: float 361 @keyword sims: A flag which if True will add the Monte Carlo simulation rotors to the structural object. There must be one model for each Monte Carlo simulation already present in the structural object. 362 @type sims: bool 363 """ 364 365 # Initialise the list structures for the rotor data. 366 axis = [] 367 span = [] 368 staggered = [] 369 pivot = [] 370 rotor_angle = [] 371 com = [] 372 label = [] 373 models = [] 374 375 # The half-rotor flag. 376 half_rotor = True 377 if representation == None: 378 half_rotor = False 379 380 # The transformation matrix (identity matrix or inversion matrix). 381 if representation == 'B': 382 T = -eye(3) 383 else: 384 T = eye(3) 385 386 # The models to loop over. 387 model_nums = [None] 388 sim_indices = [None] 389 if sims: 390 model_nums = [i+1 for i in range(cdp.sim_number)] 391 sim_indices = list(range(cdp.sim_number)) 392 393 # Loop over the models. 394 for i in range(len(model_nums)): 395 # The pivot points. 396 pivot1 = generate_pivot(order=1, sim_index=sim_indices[i], pdb_limit=True) 397 pivot2 = generate_pivot(order=2, sim_index=sim_indices[i], pdb_limit=True) 398 399 # The single rotor models. 400 if cdp.model in [MODEL_ROTOR, MODEL_FREE_ROTOR, MODEL_ISO_CONE, MODEL_ISO_CONE_FREE_ROTOR, MODEL_PSEUDO_ELLIPSE, MODEL_PSEUDO_ELLIPSE_FREE_ROTOR]: 401 # The rotor angle. 402 if cdp.model in MODEL_LIST_FREE_ROTORS: 403 rotor_angle.append(pi) 404 else: 405 if sims: 406 rotor_angle.append(cdp.cone_sigma_max_sim[sim_indices[i]]) 407 else: 408 rotor_angle.append(cdp.cone_sigma_max) 409 410 # Get the CoM of the entire molecule to use as the centre of the rotor. 411 if cdp.model in [MODEL_ROTOR, MODEL_FREE_ROTOR]: 412 com.append(pipe_centre_of_mass(verbosity=0, missing_error=False)) 413 else: 414 com.append(pivot1) 415 416 # Generate the rotor axis. 417 axes = generate_axis_system(sim_index=sim_indices[i]) 418 axis.append(axes[:, 2]) 419 420 # The size of the rotor (Angstrom), taking the cone representation into account by adding 5 Angstrom. 421 if cdp.model in [MODEL_ROTOR, MODEL_FREE_ROTOR]: 422 span.append(size) 423 else: 424 span.append(size + 5.0) 425 426 # Stagger the propeller blades. 427 if cdp.model in MODEL_LIST_FREE_ROTORS: 428 staggered.append(False) 429 else: 430 staggered.append(True) 431 432 # The pivot. 433 pivot.append(pivot1) 434 435 # The label. 436 label.append('z-ax') 437 438 # The models. 439 if sims: 440 models.append(model_nums[i]) 441 else: 442 models.append(None) 443 444 # The double rotor models. 445 elif cdp.model in [MODEL_DOUBLE_ROTOR]: 446 # Add both rotor angles (the 2nd must come first). 447 if sims: 448 rotor_angle.append(cdp.cone_sigma_max_2_sim[sim_indices[i]]) 449 rotor_angle.append(cdp.cone_sigma_max_sim[sim_indices[i]]) 450 else: 451 rotor_angle.append(cdp.cone_sigma_max_2) 452 rotor_angle.append(cdp.cone_sigma_max) 453 454 # Set the com to the pivot points. 455 com.append(pivot2) 456 com.append(pivot1) 457 458 # Generate the eigenframe of the motion. 459 frame = generate_axis_system(sim_index=sim_indices[i]) 460 461 # Add the x and y axes. 462 axis.append(frame[:, 0]) 463 axis.append(frame[:, 1]) 464 465 # The rotor size (Angstrom). 466 span.append(size) 467 span.append(size) 468 469 # Stagger the propeller blades. 470 staggered.append(True) 471 staggered.append(True) 472 473 # The pivot points. 474 pivot.append(pivot2) 475 pivot.append(pivot1) 476 477 # The labels. 478 label.append('x-ax') 479 label.append('y-ax') 480 481 # The models. 482 if sims: 483 models.append(model_nums[i]) 484 models.append(model_nums[i]) 485 else: 486 models.append(None) 487 models.append(None) 488 489 # Add each rotor to the structure as a new molecule. 490 for i in range(len(axis)): 491 rotor(structure=structure, rotor_angle=rotor_angle[i], axis=dot(T, axis[i]), axis_pt=pivot[i], label=label[i], centre=com[i], span=span[i]/1e10, blade_length=5e-10, model_num=models[i], staggered=staggered[i], half_rotor=half_rotor)

492 493

494 -def average_position(structure=None, models=None, sim=None):

495 """Shift the given structural object to the average position. 496 497 @keyword structure: The structural object to operate on. 498 @type structure: lib.structure.internal.object.Internal instance 499 @keyword models: The list of models to shift. 500 @type models: list of int 501 @keyword sim: A flag which if True will use the Monte Carlo simulation results. In this case, the model list should be set to the simulation indices plus 1 and the structural object should have one model per simulation already set up. 502 @type sim: bool 503 """ 504 505 # The selection object. 506 selection = structure.selection(atom_id=domain_moving()) 507 508 # Loop over each model. 509 for i in range(len(models)): 510 # First rotate the moving domain to the average position. 511 R = zeros((3, 3), float64) 512 if hasattr(cdp, 'ave_pos_alpha'): 513 if sim: 514 euler_to_R_zyz(cdp.ave_pos_alpha_sim[i], cdp.ave_pos_beta_sim[i], cdp.ave_pos_gamma_sim[i], R) 515 else: 516 euler_to_R_zyz(cdp.ave_pos_alpha, cdp.ave_pos_beta, cdp.ave_pos_gamma, R) 517 else: 518 if sim: 519 euler_to_R_zyz(0.0, cdp.ave_pos_beta_sim[i], cdp.ave_pos_gamma_sim[i], R) 520 else: 521 euler_to_R_zyz(0.0, cdp.ave_pos_beta, cdp.ave_pos_gamma, R) 522 origin = pipe_centre_of_mass(atom_id=domain_moving(), verbosity=0, missing_error=False) 523 structure.rotate(R=R, origin=origin, model=models[i], selection=selection) 524 525 # Then translate the moving domain. 526 if sim: 527 T = [cdp.ave_pos_x_sim[i], cdp.ave_pos_y_sim[i], cdp.ave_pos_z_sim[i]] 528 else: 529 T = [cdp.ave_pos_x, cdp.ave_pos_y, cdp.ave_pos_z] 530 structure.translate(T=T, model=models[i], selection=selection)

531 532

533 -def create_ave_pos(format='PDB', file=None, dir=None, compress_type=0, model=1, force=False):

534 """Create a PDB file of the molecule with the moving domains shifted to the average position. 535 536 @keyword format: The format for outputting the geometric representation. Currently only the 'PDB' format is supported. 537 @type format: str 538 @keyword file: The name of the file for the average molecule structure. 539 @type file: str 540 @keyword dir: The name of the directory to place the PDB file into. 541 @type dir: str 542 @keyword compress_type: The compression type. The integer values correspond to the compression type: 0, no compression; 1, Bzip2 compression; 2, Gzip compression. 543 @type compress_type: int 544 @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. 545 @type model: int 546 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten. 547 @type force: bool 548 """ 549 550 # Printout. 551 subsection(file=sys.stdout, text="Creating a PDB file with the moving domains shifted to the average position.") 552 553 # Checks. 554 if not hasattr(cdp, 'structure'): 555 warn(RelaxWarning("No structural data is present, cannot create the average position representation.")) 556 return 557 558 # Initialise. 559 titles = [] 560 sims = [] 561 file_root = [] 562 models = [] 563 structures = [] 564 565 # The real average position. 566 titles.append("real average position") 567 sims.append(False) 568 file_root.append(file) 569 models.append([None]) 570 571 # The positive MC simulation representation. 572 if hasattr(cdp, 'sim_number'): 573 titles.append("MC simulation representation") 574 sims.append(True) 575 file_root.append("%s_sim" % file) 576 models.append([i+1 for i in range(cdp.sim_number)]) 577 578 # Make a copy of the structural object (so as to preserve the original structure). 579 structures.append(deepcopy(cdp.structure)) 580 if hasattr(cdp, 'sim_number'): 581 structures.append(deepcopy(cdp.structure)) 582 583 # Delete all but the chosen model for the simulations. 584 if hasattr(cdp, 'sim_number') and len(structures[-1].structural_data) > 1: 585 # Determine the models to delete. 586 to_delete = [] 587 for model_cont in structures[-1].model_loop(): 588 if model_cont.num != model: 589 to_delete.append(model_cont.num) 590 to_delete.reverse() 591 592 # Delete them. 593 for num in to_delete: 594 structures[-1].structural_data.delete_model(model_num=num) 595 596 # Loop over each representation and add the contents. 597 for i in range(len(titles)): 598 # Printout. 599 subsubsection(file=sys.stdout, text="Creating the %s." % titles[i]) 600 601 # Loop over each model. 602 for j in range(len(models[i])): 603 # Create or set the models, if needed. 604 if models[i][j] == 1: 605 structures[i].set_model(model_new=1) 606 elif models[i][j] != None: 607 structures[i].add_model(model=models[i][j]) 608 609 # Shift to the average position. 610 average_position(structure=structures[i], models=models[i], sim=sims[i]) 611 612 # Output to PDB format. 613 if format == 'PDB': 614 pdb_file = open_write_file(file_name=file_root[i]+'.pdb', dir=dir, compress_type=compress_type, force=force) 615 structures[i].write_pdb(file=pdb_file) 616 pdb_file.close()

617 618

619 -def create_geometric_rep(format='PDB', file=None, dir=None, compress_type=0, size=30.0, inc=36, force=False):

620 """Create a PDB file containing a geometric object representing the frame order dynamics. 621 622 @keyword format: The format for outputting the geometric representation. Currently only the 'PDB' format is supported. 623 @type format: str 624 @keyword file: The name of the file of the PDB representation of the frame order dynamics to create. 625 @type file: str 626 @keyword dir: The name of the directory to place the PDB file into. 627 @type dir: str 628 @keyword compress_type: The compression type. The integer values correspond to the compression type: 0, no compression; 1, Bzip2 compression; 2, Gzip compression. 629 @type compress_type: int 630 @keyword size: The size of the geometric object in Angstroms. 631 @type size: float 632 @keyword inc: The number of increments for the filling of the cone objects. 633 @type inc: int 634 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten. 635 @type force: bool 636 """ 637 638 # Printout. 639 subsection(file=sys.stdout, text="Creating a PDB file containing a geometric object representing the frame order dynamics.") 640 641 # Checks. 642 check_parameters(escalate=2) 643 644 # Initialise. 645 titles = [] 646 structures = [] 647 representation = [] 648 sims = [] 649 file_root = [] 650 651 # Symmetry for inverted representations? 652 sym = True 653 if cdp.model in [MODEL_ROTOR, MODEL_FREE_ROTOR, MODEL_DOUBLE_ROTOR]: 654 sym = False 655 656 # The standard representation. 657 titles.append("Representation A") 658 structures.append(Internal()) 659 if sym: 660 representation.append('A') 661 file_root.append("%s_A" % file) 662 else: 663 representation.append(None) 664 file_root.append(file) 665 sims.append(False) 666 667 # The inverted representation. 668 if sym: 669 titles.append("Representation A") 670 structures.append(Internal()) 671 representation.append('B') 672 file_root.append("%s_B" % file) 673 sims.append(False) 674 675 # The standard MC simulation representation. 676 if hasattr(cdp, 'sim_number'): 677 titles.append("MC simulation representation A") 678 structures.append(Internal()) 679 if sym: 680 representation.append('A') 681 file_root.append("%s_sim_A" % file) 682 else: 683 representation.append(None) 684 file_root.append("%s_sim" % file) 685 sims.append(True) 686 687 # The inverted MC simulation representation. 688 if hasattr(cdp, 'sim_number') and sym: 689 titles.append("MC simulation representation B") 690 structures.append(Internal()) 691 representation.append('B') 692 file_root.append("%s_sim_B" % file) 693 sims.append(True) 694 695 # Loop over each structure and add the contents. 696 for i in range(len(structures)): 697 # Printout. 698 subsubsection(file=sys.stdout, text="Creating the %s." % titles[i]) 699 700 # Create a model for each Monte Carlo simulation. 701 if sims[i]: 702 for sim_i in range(cdp.sim_number): 703 structures[i].add_model(model=sim_i+1) 704 705 # Add the pivots. 706 add_pivots(structure=structures[i], sims=sims[i]) 707 708 # Add all rotor objects. 709 add_rotors(structure=structures[i], representation=representation[i], size=size, sims=sims[i]) 710 711 # Add the axis systems. 712 add_axes(structure=structures[i], representation=representation[i], size=size, sims=sims[i]) 713 714 # Add the cone objects. 715 if cdp.model not in [MODEL_ROTOR, MODEL_FREE_ROTOR, MODEL_DOUBLE_ROTOR]: 716 add_cones(structure=structures[i], representation=representation[i], size=size, inc=inc, sims=sims[i]) 717 718 # Add atoms for creating titles. 719 add_titles(structure=structures[i], representation=representation[i], displacement=size+10, sims=sims[i]) 720 721 # Create the PDB file. 722 if format == 'PDB': 723 pdb_file = open_write_file(file_root[i]+'.pdb', dir, compress_type=compress_type, force=force) 724 structures[i].write_pdb(pdb_file) 725 pdb_file.close()

726 727

728 -def frame_from_axis(axis, frame):

729 """Build a full 3D frame from the single axis. 730 731 @param axis: The Z-axis of the system. 732 @type axis: numpy rank-1, 3D array 733 @param frame: The empty frame to populate. 734 @type frame: numpy rank-2, 3D array 735 """ 736 737 # Store the Z-axis. 738 frame[:, 2] = axis 739 740 # The temporary eigenframe X-axis. 741 frame[0, 0] = 1.0 742 743 # The Y-axis (orthonormal to Z and X). 744 frame[:, 1] = cross(frame[:, 2], frame[:, 0]) 745 frame[:, 1] /= norm(frame[:, 1]) 746 747 # The orthonormal X-axis. 748 frame[:, 0] = cross(frame[:, 1], frame[:, 2]) 749 frame[:, 0] /= norm(frame[:, 0])

750 751

752 -def generate_axis_system(sim_index=None):

753 """Generate and return the full 3D axis system for the current frame order model. 754 755 @keyword sim_index: The optional MC simulation index to obtain the frame for a given simulation. 756 @type sim_index: None or int 757 @return: The full 3D axis system for the model. 758 @rtype: numpy rank-2, 3D float64 array 759 """ 760 761 # Initialise. 762 axis = zeros(3, float64) 763 frame = zeros((3, 3), float64) 764 765 # The 1st pivot point. 766 pivot = generate_pivot(order=1, sim_index=sim_index, pdb_limit=True) 767 768 # The CoM of the system. 769 com = pipe_centre_of_mass(verbosity=0, missing_error=False) 770 771 # The system for the rotor models. 772 if cdp.model in [MODEL_ROTOR, MODEL_FREE_ROTOR]: 773 # Generate the axis. 774 if sim_index == None: 775 axis = create_rotor_axis_alpha(alpha=cdp.axis_alpha, pivot=pivot, point=com) 776 else: 777 axis = create_rotor_axis_alpha(alpha=cdp.axis_alpha_sim[sim_index], pivot=pivot, point=com) 778 779 # Create a full normalised axis system. 780 frame_from_axis(axis, frame) 781 782 # The system for the isotropic cones. 783 elif cdp.model in MODEL_LIST_ISO_CONE: 784 # Generate the axis. 785 if sim_index == None: 786 axis = create_rotor_axis_spherical(theta=cdp.axis_theta, phi=cdp.axis_phi) 787 else: 788 axis = create_rotor_axis_spherical(theta=cdp.axis_theta_sim[sim_index], phi=cdp.axis_phi_sim[sim_index]) 789 790 # Create a full normalised axis system. 791 frame_from_axis(axis, frame) 792 793 # The system for the pseudo-ellipses and double rotor. 794 elif cdp.model in MODEL_LIST_PSEUDO_ELLIPSE + MODEL_LIST_DOUBLE: 795 # Generate the eigenframe of the motion. 796 if sim_index == None: 797 euler_to_R_zyz(cdp.eigen_alpha, cdp.eigen_beta, cdp.eigen_gamma, frame) 798 else: 799 euler_to_R_zyz(cdp.eigen_alpha_sim[sim_index], cdp.eigen_beta_sim[sim_index], cdp.eigen_gamma_sim[sim_index], frame) 800 801 # Unknown model. 802 else: 803 raise RelaxFault 804 805 # Return the full eigenframe. 806 return frame

807

Source Code for Module specific_analyses.frame_order.geometric