generic_fns.structure.geometric

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2013 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 # Python module imports. 23 from math import cos, pi, sin 24 from numpy import arccos, array, dot, eye, float64, transpose, zeros 25 from os import getcwd 26 from string import ascii_uppercase 27 from warnings import warn 28 29 # relax module imports. 30 from generic_fns.interatomic import interatomic_loop 31 from generic_fns.mol_res_spin import exists_mol_res_spin_data, return_spin 32 from generic_fns import pipes 33 from generic_fns.structure.internal import Internal 34 from generic_fns.structure.mass import centre_of_mass 35 from lib.structure.rotor import rotor_pdb 36 from maths_fns.rotation_matrix import two_vect_to_R 37 from relax_errors import RelaxNoPdbError, RelaxNoSequenceError, RelaxNoTensorError, RelaxNoVectorsError 38 from relax_io import get_file_path, open_write_file 39 from relax_warnings import RelaxWarning 40 from status import Status; status = Status() 41 42

43 -def angles_regular(inc=None):

44 """Determine the spherical angles for a regular sphere point distribution. 45 46 @keyword inc: The number of increments in the distribution. 47 @type inc: int 48 @return: The phi angle array and the theta angle array. 49 @rtype: array of float, array of float 50 """ 51 52 # Generate the increment values of u. 53 u = zeros(inc, float64) 54 val = 1.0 / float(inc) 55 for i in range(inc): 56 u[i] = float(i) * val 57 58 # Generate the increment values of v. 59 v = zeros(inc/2+1, float64) 60 val = 1.0 / float(inc/2) 61 for i in range(int(inc/2+1)): 62 v[i] = float(i) * val 63 64 # Generate the distribution of spherical angles theta. 65 theta = 2.0 * pi * u 66 67 # Generate the distribution of spherical angles phi (from bottom to top). 68 phi = zeros(len(v), float64) 69 for i in range(len(v)): 70 phi[len(v)-1-i] = pi * v[i] 71 72 # Return the angle arrays. 73 return phi, theta

74 75

76 -def angles_uniform(inc=None):

77 """Determine the spherical angles for a uniform sphere point distribution. 78 79 @keyword inc: The number of increments in the distribution. 80 @type inc: int 81 @return: The phi angle array and the theta angle array. 82 @rtype: array of float, array of float 83 """ 84 85 # Generate the increment values of u. 86 u = zeros(inc, float64) 87 val = 1.0 / float(inc) 88 for i in range(inc): 89 u[i] = float(i) * val 90 91 # Generate the increment values of v. 92 v = zeros(inc/2+2, float64) 93 val = 1.0 / float(inc/2) 94 for i in range(1, int(inc/2)+1): 95 v[i] = float(i-1) * val + val/2.0 96 v[-1] = 1.0 97 98 # Generate the distribution of spherical angles theta. 99 theta = 2.0 * pi * u 100 101 # Generate the distribution of spherical angles phi. 102 phi = arccos(2.0 * v - 1.0) 103 104 # Return the angle arrays. 105 return phi, theta

106 107

108 -def autoscale_tensor(method='mass'):

109 """Automatically determine an appropriate scaling factor for display of the diffusion tensor. 110 111 @keyword method: The method used to determine the scaling of the diffusion tensor object. 112 @type method: str 113 @return: The scaling factor. 114 @rtype: float 115 """ 116 117 # Centre of mass method. 118 if method == 'mass': 119 com, mass = centre_of_mass(return_mass=True) 120 scale = mass * 6.8e-11 121 return scale 122 123 # Autoscaling method. 124 warn(RelaxWarning("Autoscale method %s not implimented. Reverting to scale=1.8e-6 A.s" % method)) 125 return 1.8e-6

126 127

128 -def cone_edge(mol=None, cone=None, res_name='CON', res_num=None, chain_id='', apex=None, axis=None, R=None, scale=None, inc=None, distribution='uniform', debug=False):

129 """Add a residue to the atomic data representing a cone of the given angle. 130 131 A series of vectors totalling the number of increments and starting at the origin are equally spaced around the cone axis. The atoms representing neighbouring vectors will be directly bonded together. This will generate an object representing the outer edge of a cone. 132 133 134 @keyword mol: The molecule container. 135 @type mol: MolContainer instance 136 @keyword cone: The cone object. This should provide the limit_check() method with determines the limits of the distribution accepting two arguments, the polar angle phi and the azimuthal angle theta, and return True if the point is in the limits or False if outside. It should also provide the phi_max() method for returning the phi value for the given theta. 137 @type cone: class instance 138 @keyword res_name: The residue name. 139 @type res_name: str 140 @keyword res_num: The residue number. 141 @type res_num: int 142 @keyword chain_id: The chain identifier. 143 @type chain_id: str 144 @keyword apex: The apex of the cone. 145 @type apex: numpy array, len 3 146 @keyword axis: The central axis of the cone. If supplied, then this arg will be used to construct the rotation matrix. 147 @type axis: numpy array, len 3 148 @keyword R: A 3x3 rotation matrix. If the axis arg supplied, then this matrix will be ignored. 149 @type R: 3x3 numpy array 150 @keyword scale: The scaling factor to stretch all points by. 151 @type scale: float 152 @keyword inc: The number of increments or number of vectors used to generate the outer edge of the cone. 153 @type inc: int 154 @keyword distribution: The type of point distribution to use. This can be 'uniform' or 'regular'. 155 @type distribution: str 156 """ 157 158 # The atom numbers (and indices). 159 atom_num = mol.atom_num[-1]+1 160 161 # Add an atom for the cone apex. 162 mol.atom_add(pdb_record='HETATM', atom_num=atom_num, atom_name='APX', res_name=res_name, res_num=res_num, pos=apex, segment_id=None, element='H') 163 origin_atom = atom_num 164 165 # Get the polar and azimuthal angles for the distribution. 166 if distribution == 'uniform': 167 phi, theta = angles_uniform(inc) 168 else: 169 phi, theta = angles_regular(inc) 170 171 # Initialise the rotation matrix. 172 if R == None: 173 R = eye(3) 174 175 # Get the rotation matrix. 176 if axis != None: 177 two_vect_to_R(array([0, 0, 1], float64), axis, R) 178 179 # Determine the maximum phi values and the indices of the point just above the edge. 180 phi_max = zeros(len(theta), float64) 181 edge_index = zeros(len(theta), int) 182 for i in range(len(theta)): 183 # Get the polar angle for the longitude edge atom. 184 phi_max[i] = cone.phi_max(theta[i]) 185 186 # The index. 187 for j in range(len(phi)): 188 if phi[j] <= phi_max[i]: 189 edge_index[i] = j 190 break 191 192 # Reverse edge index. 193 edge_index_rev = len(phi) - 1 - edge_index 194 195 # Move around the azimuth. 196 atom_num = atom_num + 1 197 for i in range(len(theta)): 198 # The vector in the unrotated frame. 199 x = cos(theta[i]) * sin(phi_max[i]) 200 y = sin(theta[i])* sin(phi_max[i]) 201 z = cos(phi_max[i]) 202 vector = array([x, y, z], float64) 203 204 # Rotate the vector. 205 vector = dot(R, vector) 206 207 # The atom id. 208 atom_id = 'T' + repr(i) 209 210 # The atom position. 211 pos = apex+vector*scale 212 213 # Debugging. 214 if debug: 215 print("i: %s; theta: %s" % (i, theta[i])) 216 print("%sAdding atom %s." % (" "*4, get_proton_name(atom_num))) 217 218 # Add the vector as a H atom of the cone residue. 219 mol.atom_add(pdb_record='HETATM', atom_num=atom_num, atom_name=get_proton_name(atom_num), res_name=res_name, res_num=res_num, pos=pos, segment_id=None, element='H') 220 221 # Join the longitude atom to the cone apex. 222 mol.atom_connect(index1=origin_atom-1, index2=atom_num-1) 223 224 # Increment the atom number. 225 atom_num = atom_num + 1 226 227 # Add latitude atoms for a smoother cone edge and better stitching to the cap. 228 for j in range(len(phi)): 229 # The index for the direction top to bottom! 230 k = len(phi) - 1 - j 231 232 # Debugging. 233 if debug: 234 print("%sj: %s; phi: %-20s; k: %s; phi: %-20s; phi_max: %-20s" % (" "*4, j, phi[j], k, phi[k], phi_max[i])) 235 236 # No edge. 237 skip = True 238 239 # Forward edge (skip when the latitude is phi max). 240 fwd_index = i+1 241 if i == len(theta)-1: 242 fwd_index = 0 243 if j >= edge_index[i] and j < edge_index[fwd_index] and not abs(phi_max[fwd_index] - phi[j]) < 1e-6: 244 # Debugging. 245 if debug: 246 print("%sForward edge." % (" "*8)) 247 248 # Edge found. 249 skip = False 250 251 # Find the theta value for this polar angle phi. 252 phi_val = phi[j] 253 if fwd_index == 0: 254 theta_max = theta[fwd_index] + 2*pi 255 else: 256 theta_max = theta[fwd_index] 257 theta_max = cone.theta_max(phi_val, theta_min=theta[i], theta_max=theta_max-1e-7) 258 259 # Back edge (skip when the latitude is phi max). 260 rev_index = i-1 261 if i == 0: 262 rev_index = len(theta)-1 263 if i < len(theta)-1 and j > edge_index_rev[i] and j <= edge_index_rev[i+1] and not abs(phi_max[fwd_index] - phi[k]) < 1e-6: 264 # Debugging. 265 if debug: 266 print("%sBack edge." % (" "*8)) 267 268 # Edge found. 269 skip = False 270 271 # Find the theta value for this polar angle phi. 272 phi_val = phi[k] 273 theta_max = cone.theta_max(phi_val, theta_min=theta[i-1], theta_max=theta[i]) 274 275 # Skipping. 276 if skip: 277 continue 278 279 # Debugging. 280 if debug: 281 print("%sAdding atom %s." % (" "*8, get_proton_name(atom_num))) 282 283 # The coordinates. 284 x = cos(theta_max) * sin(phi_val) 285 y = sin(theta_max) * sin(phi_val) 286 z = cos(phi_val) 287 pos = array([x, y, z], float64) * scale 288 289 # Rotate and shift. 290 pos = apex + dot(R, pos) 291 292 # Add the atom. 293 mol.atom_add(pdb_record='HETATM', atom_num=atom_num, atom_name=get_proton_name(atom_num), res_name=res_name, chain_id=chain_id, res_num=res_num, pos=pos, segment_id=None, element='H') 294 295 # Increment the atom number. 296 atom_num = atom_num + 1 297 298 # Debugging. 299 if debug: 300 print("\nBuilding the edge.") 301 302 # Build the cone edge. 303 for i in range(origin_atom, atom_num-2): 304 # Debugging. 305 if debug: 306 print("%sCone edge, connecting %s to %s" % (" "*4, get_proton_name(i), get_proton_name(i+1))) 307 308 # Connect. 309 mol.atom_connect(index1=i, index2=i+1) 310 311 # Connect the last radial array to the first (to zip up the circle). 312 mol.atom_connect(index1=atom_num-2, index2=origin_atom)

313 314

315 -def create_cone_pdb(mol=None, cone=None, start_res=1, apex=None, axis=None, R=None, inc=None, scale=30.0, distribution='regular', file=None, dir=None, force=False, axis_flag=True):

316 """Create a PDB representation of the given cone object. 317 318 @keyword mol: The molecule container. 319 @type mol: MolContainer instance 320 @keyword cone: The cone object. This should provide the limit_check() method with determines the limits of the distribution accepting two arguments, the polar angle phi and the azimuthal angle theta, and return True if the point is in the limits or False if outside. It should also provide the theta_max() method for returning the theta value for the given phi, the phi_max() method for returning the phi value for the given theta. 321 @type cone: class instance 322 @keyword start_res: The starting residue number. 323 @type start_res: str 324 @keyword apex: The apex of the cone. 325 @type apex: rank-1, 3D numpy array 326 @keyword axis: The central axis of the cone. If not supplied, the z-axis will be used. 327 @type axis: rank-1, 3D numpy array 328 @keyword R: The rotation matrix. 329 @type R: rank-2, 3D numpy array 330 @keyword inc: The increment number used to determine the number of latitude and longitude lines. 331 @type inc: int 332 @keyword scale: The scaling factor to stretch the unit cone by. 333 @type scale: float 334 @keyword distribution: The type of point distribution to use. This can be 'uniform' or 'regular'. 335 @type distribution: str 336 @keyword file: The name of the PDB file to create. 337 @type file: str 338 @keyword dir: The name of the directory to place the PDB file into. 339 @type dir: str 340 @keyword force: Flag which if set to True will overwrite any pre-existing file. 341 @type force: bool 342 @keyword axis_flag: A flag which if True will create the cone's axis. 343 @type axis_flag: bool 344 """ 345 346 # The cone axis default of the z-axis. 347 if not axis: 348 axis = array([0, 0, 1], float64) 349 350 # No rotation. 351 if R == None: 352 R = eye(3) 353 354 # No molecule supplied. 355 if mol == None: 356 # Create the structural object. 357 structure = Internal() 358 359 # Add a molecule. 360 structure.add_molecule(name='cone') 361 362 # Alias the single molecule from the single model. 363 mol = structure.structural_data[0].mol[0] 364 365 # The first atom number. 366 start_atom = 1 367 if hasattr(mol, 'atom_num'): 368 start_atom = mol.atom_num[-1]+1 369 370 # The axis. 371 if axis_flag: 372 # Add the apex. 373 mol.atom_add(pdb_record='HETATM', atom_num=start_atom, atom_name='R', res_name='APX', res_num=start_res, pos=apex, element='C') 374 375 # Generate the axis vectors. 376 print("\nGenerating the axis vectors.") 377 res_num = generate_vector_residues(mol=mol, vector=dot(R, axis), atom_name='Axis', res_name_vect='AXE', res_num=start_res+1, origin=apex, scale=scale) 378 379 # Generate the cone outer edge. 380 print("\nGenerating the cone outer edge.") 381 edge_start_atom = mol.atom_num[-1]+1 382 cone_edge(mol=mol, cone=cone, res_name='EDG', res_num=start_res+2, apex=apex, R=R, scale=scale, inc=inc, distribution=distribution) 383 384 # Generate the cone cap, and stitch it to the cone edge. 385 print("\nGenerating the cone cap.") 386 cone_start_atom = mol.atom_num[-1]+1 387 generate_vector_dist(mol=mol, res_name='CON', res_num=start_res+3, centre=apex, R=R, limit_check=cone.limit_check, scale=scale, inc=inc, distribution=distribution) 388 stitch_cone_to_edge(mol=mol, cone=cone, dome_start=cone_start_atom, edge_start=edge_start_atom+1, scale=scale, inc=inc, distribution=distribution) 389 390 # Create the PDB file. 391 if file != None: 392 print("\nGenerating the PDB file.") 393 pdb_file = open_write_file(file_name=file, dir=dir, force=force) 394 structure.write_pdb(pdb_file) 395 pdb_file.close() 396 397 # Add the file to the results file list. 398 if not hasattr(cdp, 'result_files'): 399 cdp.result_files = [] 400 cdp.result_files.append(['cone_pdb', 'Cone PDB', get_file_path(file, dir)]) 401 status.observers.result_file.notify()

402 403

404 -def create_diff_tensor_pdb(scale=1.8e-6, file=None, dir=None, force=False):

405 """Create the PDB representation of the diffusion tensor. 406 407 @keyword scale: The scaling factor for the diffusion tensor. 408 @type scale: float 409 @keyword file: The name of the PDB file to create. 410 @type file: str 411 @keyword dir: The name of the directory to place the PDB file into. 412 @type dir: str 413 @keyword force: Flag which if set to True will overwrite any pre-existing file. 414 @type force: bool 415 """ 416 417 # Arguments. 418 if scale == 'mass': 419 scale = autoscale_tensor(scale) 420 421 # Test if the current data pipe exists. 422 pipes.test() 423 424 # Create an array of data pipes to loop over (hybrid support). 425 if cdp.pipe_type == 'hybrid': 426 pipe_list = cdp.hybrid_pipes 427 else: 428 pipe_list = [pipes.cdp_name()] 429 430 # Create the structural object. 431 structure = Internal() 432 433 # Add a structure. 434 structure.add_molecule(name='diff_tensor') 435 436 # Alias the single molecule from the single model. 437 mol = structure.get_molecule('diff_tensor') 438 439 # Loop over the pipes. 440 for pipe_index in range(len(pipe_list)): 441 # Get the pipe container. 442 pipe = pipes.get_pipe(pipe_list[pipe_index]) 443 444 445 # Tests. 446 ######## 447 448 # Test if the diffusion tensor data is loaded. 449 if not hasattr(pipe, 'diff_tensor'): 450 raise RelaxNoTensorError('diffusion') 451 452 # Test if a structure has been loaded. 453 if not hasattr(cdp, 'structure'): 454 raise RelaxNoPdbError 455 456 457 # Initialise. 458 ############# 459 460 # Initialise the residue number. 461 res_num = 1 462 463 # The chain identifier. 464 chain_id = ascii_uppercase[pipe_index] 465 466 # Atom ID extension (allow for multiple chains for hybrid data pipes). 467 atom_id_ext = '_' + chain_id 468 469 # Print out. 470 print("\nChain " + chain_id + "\n") 471 472 473 # Centre of mass. 474 ################# 475 476 # Calculate the centre of mass. 477 CoM = centre_of_mass() 478 479 # Add the central atom. 480 mol.atom_add(pdb_record='HETATM', atom_num=1, atom_name='R'+atom_id_ext, res_name='COM', chain_id=chain_id, res_num=res_num, pos=CoM, segment_id=None, element='C') 481 482 # Increment the residue number. 483 res_num = res_num + 1 484 485 486 # Vector distribution. 487 ###################### 488 489 # Print out. 490 print("\nGenerating the geometric object.") 491 492 # Swap the x and z axes for the oblate spheroid (the vector distributions are centred on the z-axis). 493 if pipe.diff_tensor.type == 'spheroid' and pipe.diff_tensor.spheroid_type == 'oblate': 494 # 90 deg rotation about the diffusion frame y-axis. 495 y_rot = array([[ 0, 0, 1], 496 [ 0, 1, 0], 497 [ -1, 0, 0]], float64) 498 499 # Rotate the tensor and rotation matrix. 500 rotation = dot(transpose(pipe.diff_tensor.rotation), y_rot) 501 tensor = dot(y_rot, dot(pipe.diff_tensor.tensor_diag, transpose(y_rot))) 502 tensor = dot(rotation, dot(tensor, transpose(rotation))) 503 504 # All other systems stay as normal. 505 else: 506 tensor = pipe.diff_tensor.tensor 507 rotation = pipe.diff_tensor.rotation 508 509 # The distribution. 510 generate_vector_dist(mol=mol, res_name='TNS', res_num=res_num, chain_id=chain_id, centre=CoM, R=rotation, warp=tensor, scale=scale, inc=20) 511 512 # Increment the residue number. 513 res_num = res_num + 1 514 515 516 # Axes of the tensor. 517 ##################### 518 519 # Create the unique axis of the spheroid. 520 if pipe.diff_tensor.type == 'spheroid': 521 # Print out. 522 print("\nGenerating the unique axis of the diffusion tensor.") 523 print(" Scaling factor: " + repr(scale)) 524 525 # Simulations. 526 if hasattr(pipe.diff_tensor, 'tm_sim'): 527 sim_vectors = [] 528 for i in range(len(pipe.diff_tensor.tm_sim)): 529 sim_vectors.append(pipe.diff_tensor.Dpar_sim[i] * pipe.diff_tensor.Dpar_unit_sim[i]) 530 else: 531 sim_vectors = None 532 533 # Generate the axes representation. 534 res_num = generate_vector_residues(mol=mol, vector=pipe.diff_tensor.Dpar*pipe.diff_tensor.Dpar_unit, atom_name='Dpar', res_name_vect='AXS', sim_vectors=sim_vectors, chain_id=chain_id, res_num=res_num, origin=CoM, scale=scale, neg=True) 535 536 537 # Create the three axes of the ellipsoid. 538 if pipe.diff_tensor.type == 'ellipsoid': 539 # Print out. 540 print("Generating the three axes of the ellipsoid.") 541 print(" Scaling factor: " + repr(scale)) 542 543 # Simulations. 544 if hasattr(pipe.diff_tensor, 'tm_sim'): 545 sim_Dx_vectors = [] 546 sim_Dy_vectors = [] 547 sim_Dz_vectors = [] 548 for i in range(len(pipe.diff_tensor.tm_sim)): 549 sim_Dx_vectors.append(pipe.diff_tensor.Dx_sim[i] * pipe.diff_tensor.Dx_unit_sim[i]) 550 sim_Dy_vectors.append(pipe.diff_tensor.Dy_sim[i] * pipe.diff_tensor.Dy_unit_sim[i]) 551 sim_Dz_vectors.append(pipe.diff_tensor.Dz_sim[i] * pipe.diff_tensor.Dz_unit_sim[i]) 552 else: 553 sim_Dx_vectors = None 554 sim_Dy_vectors = None 555 sim_Dz_vectors = None 556 557 # Generate the axes representation. 558 res_num = generate_vector_residues(mol=mol, vector=pipe.diff_tensor.Dx*pipe.diff_tensor.Dx_unit, atom_name='Dx', res_name_vect='AXS', sim_vectors=sim_Dx_vectors, chain_id=chain_id, res_num=res_num, origin=CoM, scale=scale, neg=True) 559 res_num = generate_vector_residues(mol=mol, vector=pipe.diff_tensor.Dy*pipe.diff_tensor.Dy_unit, atom_name='Dy', res_name_vect='AXS', sim_vectors=sim_Dy_vectors, chain_id=chain_id, res_num=res_num, origin=CoM, scale=scale, neg=True) 560 res_num = generate_vector_residues(mol=mol, vector=pipe.diff_tensor.Dz*pipe.diff_tensor.Dz_unit, atom_name='Dz', res_name_vect='AXS', sim_vectors=sim_Dz_vectors, chain_id=chain_id, res_num=res_num, origin=CoM, scale=scale, neg=True) 561 562 563 # Create the PDB file. 564 ###################### 565 566 # Print out. 567 print("\nGenerating the PDB file.") 568 569 # Open the PDB file for writing. 570 tensor_pdb_file = open_write_file(file, dir, force=force) 571 572 # Write the data. 573 structure.write_pdb(tensor_pdb_file) 574 575 # Close the file. 576 tensor_pdb_file.close() 577 578 # Add the file to the results file list. 579 if not hasattr(cdp, 'result_files'): 580 cdp.result_files = [] 581 if dir == None: 582 dir = getcwd() 583 cdp.result_files.append(['diff_tensor_pdb', 'Diffusion tensor PDB', get_file_path(file, dir)]) 584 status.observers.result_file.notify()

585 586

587 -def create_rotor_pdb(file=None, dir=None, rotor_angle=None, axis=None, axis_pt=True, centre=None, span=2e-9, blade_length=5e-10, force=False, staggered=False):

588 """Create a PDB representation of a rotor motional model. 589 590 @keyword file: The name of the PDB file to create. 591 @type file: str 592 @keyword dir: The name of the directory to place the PDB file into. 593 @type dir: str 594 @keyword rotor_angle: The angle of the rotor motion in degrees. 595 @type rotor_angle: float 596 @keyword axis: The vector defining the rotor axis. 597 @type axis: numpy rank-1, 3D array 598 @keyword axis_pt: A point lying anywhere on the rotor axis. This is used to define the position of the axis in 3D space. 599 @type axis_pt: numpy rank-1, 3D array 600 @keyword centre: The central point of the representation. If this point is not on the rotor axis, then the closest point on the axis will be used for the centre. 601 @type centre: numpy rank-1, 3D array 602 @keyword span: The distance from the central point to the rotor blades (meters). 603 @type span: float 604 @keyword blade_length: The length of the representative rotor blades. 605 @type blade_length: float 606 @keyword force: A flag which if set will overwrite any pre-existing file. 607 @type force: bool 608 @keyword staggered: A flag which if True will cause the rotor blades to be staggered. This is used to avoid blade overlap. 609 @type staggered: bool 610 """ 611 612 # Test if the current pipe exists. 613 pipes.test() 614 615 # Convert the angle to radians. 616 rotor_angle = rotor_angle / 360.0 * 2.0 * pi 617 618 # Create the structural object. 619 structure = Internal() 620 621 # Generate the rotor object. 622 rotor_pdb(structure=structure, rotor_angle=rotor_angle, axis=axis, axis_pt=axis_pt, centre=centre, span=span, blade_length=blade_length, staggered=staggered) 623 624 # Print out. 625 print("\nGenerating the PDB file.") 626 627 # Open the PDB file for writing. 628 tensor_pdb_file = open_write_file(file, dir, force=force) 629 630 # Write the data. 631 structure.write_pdb(tensor_pdb_file) 632 633 # Close the file. 634 tensor_pdb_file.close() 635 636 # Add the file to the results file list. 637 if not hasattr(cdp, 'result_files'): 638 cdp.result_files = [] 639 if dir == None: 640 dir = getcwd() 641 cdp.result_files.append(['diff_tensor_pdb', 'Diffusion tensor PDB', get_file_path(file, dir)]) 642 status.observers.result_file.notify()

643 644

645 -def create_vector_dist(length=None, symmetry=True, file=None, dir=None, force=False):

646 """Create a PDB representation of the vector distribution. 647 648 @keyword length: The length to set the vectors to in the PDB file. 649 @type length: float 650 @keyword symmetry: The symmetry flag which if set will create a second PDB chain 'B' which is the same as chain 'A' but with the vectors reversed. 651 @type symmetry: bool 652 @keyword file: The name of the PDB file to create. 653 @type file: str 654 @keyword dir: The name of the directory to place the PDB file into. 655 @type dir: str 656 @keyword force: Flag which if set will overwrite any pre-existing file. 657 @type force: bool 658 """ 659 660 # Test if the current pipe exists. 661 pipes.test() 662 663 # Test if a structure has been loaded. 664 if not hasattr(cdp, 'structure') or not cdp.structure.num_models() > 0: 665 raise RelaxNoPdbError 666 667 # Test if sequence data is loaded. 668 if not exists_mol_res_spin_data(): 669 raise RelaxNoSequenceError 670 671 # Test if unit vectors exist. 672 vectors = False 673 for interatom in interatomic_loop(): 674 if hasattr(interatom, 'vector'): 675 vectors = True 676 break 677 if not vectors: 678 raise RelaxNoVectorsError 679 680 681 # Initialise. 682 ############# 683 684 # Create the structural object. 685 structure = Internal() 686 687 # Add a structure. 688 structure.add_molecule(name='vector_dist') 689 690 # Alias the single molecule from the single model. 691 mol = structure.structural_data[0].mol[0] 692 693 # Initialise the residue and atom numbers. 694 res_num = 1 695 atom_num = 1 696 697 698 # Centre of mass. 699 ################# 700 701 # Calculate the centre of mass. 702 R = centre_of_mass() 703 704 # Increment the residue number. 705 res_num = res_num + 1 706 707 708 # The vectors. 709 ############## 710 711 # Loop over the interatomic data containers. 712 for interatom in interatomic_loop(): 713 # Get the spins. 714 spin1 = return_spin(interatom.spin_id1) 715 spin2 = return_spin(interatom.spin_id2) 716 717 # Skip deselected spin systems. 718 if not spin1.select or not spin2.select: 719 continue 720 721 # Skip containers missing vectors. 722 if not hasattr(interatom, 'vector'): 723 continue 724 725 # Scale the vector. 726 vector = interatom.vector * length * 1e10 727 728 # Add the first spin as the central atom. 729 mol.atom_add(pdb_record='ATOM', atom_num=atom_num, atom_name=spin1.name, res_name=spin1._res_name, chain_id='A', res_num=spin1._res_num, pos=R, segment_id=None, element=spin1.element) 730 731 # Add the second spin as the end atom. 732 mol.atom_add(pdb_record='ATOM', atom_num=atom_num+1, atom_name=spin2.name, res_name=spin2._res_name, chain_id='A', res_num=spin2._res_num, pos=R+vector, segment_id=None, element=spin2.element) 733 734 # Connect the two atoms. 735 mol.atom_connect(index1=atom_num-1, index2=atom_num) 736 737 # Increment the atom number. 738 atom_num = atom_num + 2 739 740 # Symmetry chain. 741 if symmetry: 742 # Loop over the interatomic data containers. 743 for interatom in interatomic_loop(): 744 # Get the spins. 745 spin1 = return_spin(interatom.spin_id1) 746 spin2 = return_spin(interatom.spin_id2) 747 748 # Skip deselected spin systems. 749 if not spin1.select or not spin2.select: 750 continue 751 752 # Skip containers missing vectors. 753 if not hasattr(interatom, 'vector'): 754 continue 755 756 # Scale the vector. 757 vector = interatom.vector * length * 1e10 758 759 # Add the first spin as the central atom. 760 mol.atom_add(pdb_record='ATOM', atom_num=atom_num, atom_name=spin1.name, res_name=spin1._res_name, chain_id='B', res_num=spin1._res_num, pos=R, segment_id=None, element=spin1.element) 761 762 # Add the second spin as the end atom. 763 mol.atom_add(pdb_record='ATOM', atom_num=atom_num+1, atom_name=spin2.name, res_name=spin2._res_name, chain_id='B', res_num=spin2._res_num, pos=R-vector, segment_id=None, element=spin2.element) 764 765 # Connect the two atoms. 766 mol.atom_connect(index1=atom_num-1, index2=atom_num) 767 768 # Increment the atom number. 769 atom_num = atom_num + 2 770 771 772 # Create the PDB file. 773 ###################### 774 775 # Print out. 776 print("\nGenerating the PDB file.") 777 778 # Open the PDB file for writing. 779 tensor_pdb_file = open_write_file(file, dir, force=force) 780 781 # Write the data. 782 structure.write_pdb(tensor_pdb_file) 783 784 # Close the file. 785 tensor_pdb_file.close() 786 787 # Add the file to the results file list. 788 if not hasattr(cdp, 'result_files'): 789 cdp.result_files = [] 790 if dir == None: 791 dir = getcwd() 792 cdp.result_files.append(['vector_dist_pdb', 'Vector distribution PDB', get_file_path(file, dir)]) 793 status.observers.result_file.notify()

794 795

796 -def generate_vector_dist(mol=None, res_name=None, res_num=None, chain_id='', centre=zeros(3, float64), R=eye(3), warp=eye(3), limit_check=None, scale=1.0, inc=20, distribution='uniform', debug=False):

797 """Generate a uniformly distributed distribution of atoms on a warped sphere. 798 799 The vectors from the function vect_dist_spherical_angles() are used to generate the distribution. These vectors are rotated to the desired frame using the rotation matrix 'R', then each compressed or stretched by the dot product with the 'warp' matrix. Each vector is centred and at the head of the vector, a proton is placed. 800 801 802 @keyword mol: The molecule container. 803 @type mol: MolContainer instance 804 @keyword res_name: The residue name. 805 @type res_name: str 806 @keyword res_num: The residue number. 807 @type res_num: int 808 @keyword chain_id: The chain identifier. 809 @type chain_id: str 810 @keyword centre: The centre of the distribution. 811 @type centre: numpy array, len 3 812 @keyword R: The optional 3x3 rotation matrix. 813 @type R: 3x3 numpy array 814 @keyword warp: The optional 3x3 warping matrix. 815 @type warp: 3x3 numpy array 816 @keyword limit_check: A function with determines the limits of the distribution. It should accept two arguments, the polar angle phi and the azimuthal angle theta, and return True if the point is in the limits or False if outside. 817 @type limit_check: function 818 @keyword scale: The scaling factor to stretch all rotated and warped vectors by. 819 @type scale: float 820 @keyword inc: The number of increments or number of vectors used to generate the outer edge of the cone. 821 @type inc: int 822 @keyword distribution: The type of point distribution to use. This can be 'uniform' or 'regular'. 823 @type distribution: str 824 """ 825 826 # Initial atom number. 827 if len(mol.atom_num) == 0: 828 origin_num = 1 829 else: 830 origin_num = mol.atom_num[-1]+1 831 atom_num = origin_num 832 833 # Get the uniform vector distribution. 834 print(" Creating the uniform vector distribution.") 835 vectors = vect_dist_spherical_angles(inc=inc, distribution=distribution) 836 837 # Get the polar and azimuthal angles for the distribution. 838 if distribution == 'uniform': 839 phi, theta = angles_uniform(inc) 840 else: 841 phi, theta = angles_regular(inc) 842 843 # Init the arrays for stitching together. 844 edge = zeros(len(theta)) 845 edge_index = zeros(len(theta), int) 846 edge_phi = zeros(len(theta), float64) 847 edge_atom = zeros(len(theta)) 848 849 # Loop over the radial array of vectors (change in longitude). 850 for i in range(len(theta)): 851 # Debugging. 852 if debug: 853 print("i: %s; theta: %s" % (i, theta[i])) 854 855 # Loop over the vectors of the radial array (change in latitude). 856 for j in range(len(phi)): 857 # Debugging. 858 if debug: 859 print("%sj: %s; phi: %s" % (" "*4, j, phi[j])) 860 861 # Skip the vector if the point is outside of the limits. 862 if limit_check and not limit_check(phi[j], theta[i]): 863 if debug: 864 print("%sOut of cone." % (" "*8)) 865 continue 866 867 # Update the edge for this longitude. 868 if not edge[i]: 869 edge[i] = 1 870 edge_index[i] = j 871 edge_phi[i] = phi[j] 872 edge_atom[i] = atom_num 873 874 # Debugging. 875 if debug: 876 print("%sEdge detected." % (" "*8)) 877 print("%sEdge index: %s" % (" "*8, edge_index[i])) 878 print("%sEdge phi pos: %s" % (" "*8, edge_phi[i])) 879 print("%sEdge atom: %s" % (" "*8, edge_atom[i])) 880 881 # Rotate the vector into the diffusion frame. 882 vector = dot(R, vectors[i + j*len(theta)]) 883 884 # Set the length of the vector to its diffusion rate within the diffusion tensor geometric object. 885 vector = dot(warp, vector) 886 887 # Scale the vector. 888 vector = vector * scale 889 890 # Position relative to the centre of mass. 891 pos = centre + vector 892 893 # Debugging. 894 if debug: 895 print("%sAdding atom %s." % (" "*8, get_proton_name(atom_num))) 896 897 # Add the vector as a H atom of the TNS residue. 898 mol.atom_add(pdb_record='HETATM', atom_num=atom_num, atom_name=get_proton_name(atom_num), res_name=res_name, chain_id=chain_id, res_num=res_num, pos=pos, segment_id=None, element='H') 899 900 # Connect to the previous atom to generate the longitudinal lines (except for the first point). 901 if j > edge_index[i]: 902 # Debugging. 903 if debug: 904 print("%sLongitude line, connecting %s to %s" % (" "*8, get_proton_name(atom_num), get_proton_name(atom_num-1))) 905 906 mol.atom_connect(index1=atom_num-1, index2=atom_num-2) 907 908 # Connect across the radial arrays (to generate the latitudinal lines). 909 if i != 0 and j >= edge_index[i-1]: 910 # The number of atoms back to the previous longitude. 911 num = len(phi) - edge_index[i] 912 913 # Debugging. 914 if debug: 915 print("%sLatitude line, connecting %s to %s" % (" "*8, get_proton_name(atom_num), get_proton_name(atom_num-num))) 916 917 mol.atom_connect(index1=atom_num-1, index2=atom_num-1-num) 918 919 # Connect the last radial array to the first (to zip up the geometric object and close the latitudinal lines). 920 if i == len(theta)-1 and j >= edge_index[0]: 921 # The number of atom in the first longitude line. 922 num = origin_num + j - edge_index[0] 923 924 # Debugging. 925 if debug: 926 print("%sZipping up, connecting %s to %s" % (" "*8, get_proton_name(atom_num), get_proton_name(num))) 927 928 mol.atom_connect(index1=atom_num-1, index2=num-1) 929 930 # Increment the atom number. 931 atom_num = atom_num + 1

932 933

934 -def generate_vector_residues(mol=None, vector=None, atom_name=None, res_name_vect='AXS', sim_vectors=None, res_name_sim='SIM', chain_id='', res_num=None, origin=None, scale=1.0, label_placement=1.1, neg=False):

935 """Generate residue representations for the vector and the MC simulationed vectors. 936 937 This is used to create a PDB representation of any vector, including its Monte Carlo simulations. 938 939 @keyword mol: The molecule container. 940 @type mol: MolContainer instance 941 @keyword vector: The vector to be represented in the PDB. 942 @type vector: numpy array, len 3 943 @keyword atom_name: The atom name used to label the atom representing the head of the vector. 944 @type atom_name: str 945 @keyword res_name_vect: The 3 letter PDB residue code used to represent the vector. 946 @type res_name_vect: str 947 @keyword sim_vectors: The optional Monte Carlo simulation vectors to be represented in the PDB. 948 @type sim_vectors: list of numpy array, each len 3 949 @keyword res_name_sim: The 3 letter PDB residue code used to represent the Monte Carlo simulation vectors. 950 @type res_name_sim: str 951 @keyword chain_id: The chain identification code. 952 @type chain_id: str 953 @keyword res_num: The residue number. 954 @type res_num: int 955 @keyword origin: The origin for the axis. 956 @type origin: numpy array, len 3 957 @keyword scale: The scaling factor to stretch the vectors by. 958 @type scale: float 959 @keyword label_placement: A scaling factor to multiply the pre-scaled vector by. This is used to place the vector labels a little further out from the vector itself. 960 @type label_placement: float 961 @keyword neg: If True, then the negative vector positioned at the origin will also be included. 962 @type neg: bool 963 @return: The new residue number. 964 @rtype: int 965 """ 966 967 # The atom numbers (and indices). 968 origin_num = mol.atom_num[-1]+1 969 atom_num = mol.atom_num[-1]+2 970 atom_neg_num = mol.atom_num[-1]+3 971 972 # The origin atom. 973 mol.atom_add(pdb_record='HETATM', atom_num=origin_num, atom_name='R', res_name=res_name_vect, chain_id=chain_id, res_num=res_num, pos=origin, segment_id=None, element='C') 974 975 # Create the PDB residue representing the vector. 976 mol.atom_add(pdb_record='HETATM', atom_num=atom_num, atom_name=atom_name, res_name=res_name_vect, chain_id=chain_id, res_num=res_num, pos=origin+vector*scale, segment_id=None, element='C') 977 mol.atom_connect(index1=atom_num-1, index2=origin_num-1) 978 num = 1 979 if neg: 980 mol.atom_add(pdb_record='HETATM', atom_num=atom_neg_num, atom_name=atom_name, res_name=res_name_vect, chain_id=chain_id, res_num=res_num, pos=origin-vector*scale, segment_id=None, element='C') 981 mol.atom_connect(index1=atom_neg_num-1, index2=origin_num-1) 982 num = 2 983 984 # Add another atom to allow the axis labels to be shifted just outside of the vector itself. 985 mol.atom_add(pdb_record='HETATM', atom_num=atom_num+num, atom_name=atom_name, res_name=res_name_vect, chain_id=chain_id, res_num=res_num, pos=origin+label_placement*vector*scale, segment_id=None, element='N') 986 if neg: 987 mol.atom_add(pdb_record='HETATM', atom_num=atom_neg_num+num, atom_name=atom_name, res_name=res_name_vect, chain_id=chain_id, res_num=res_num, pos=origin-label_placement*vector*scale, segment_id=None, element='N') 988 989 # Print out. 990 print(" " + atom_name + " vector (scaled + shifted to origin): " + repr(origin+vector*scale)) 991 print(" Creating the MC simulation vectors.") 992 993 # Monte Carlo simulations. 994 if sim_vectors != None: 995 for i in range(len(sim_vectors)): 996 # Increment the residue number, so each simulation is a new residue. 997 res_num = res_num + 1 998 999 # The atom numbers (and indices). 1000 atom_num = mol.atom_num[-1]+1 1001 atom_neg_num = mol.atom_num[-1]+2 1002 1003 # Create the PDB residue representing the vector. 1004 mol.atom_add(pdb_record='HETATM', atom_num=atom_num, atom_name=atom_name, res_name=res_name_sim, chain_id=chain_id, res_num=res_num, pos=origin+sim_vectors[i]*scale, segment_id=None, element='C') 1005 mol.atom_connect(index1=atom_num-1, index2=origin_num-1) 1006 if neg: 1007 mol.atom_add(pdb_record='HETATM', atom_num=atom_num+1, atom_name=atom_name, res_name=res_name_sim, chain_id=chain_id, res_num=res_num, pos=origin-sim_vectors[i]*scale, segment_id=None, element='C') 1008 mol.atom_connect(index1=atom_neg_num-1, index2=origin_num-1) 1009 1010 # Return the new residue number. 1011 return res_num

1012 1013

1014 -def get_proton_name(atom_num):

1015 """Return a valid PDB atom name of <4 characters. 1016 1017 @param atom_num: The number of the atom. 1018 @type atom_num: int 1019 @return: The atom name to use in the PDB. 1020 @rtype: str 1021 """ 1022 1023 # Init the proton first letters and the atom number folding limits. 1024 names = ['H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q'] 1025 lims = [0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000] 1026 1027 # Loop over the proton names. 1028 for i in range(len(names)): 1029 # In the bounds. 1030 if atom_num >= lims[i] and atom_num < lims[i+1]: 1031 return names[i] + repr(atom_num - lims[i])

1032 1033

1034 -def stitch_cone_to_edge(mol=None, cone=None, chain_id='', dome_start=None, edge_start=None, scale=1.0, inc=None, distribution='uniform', debug=False):

1035 """Function for stitching the cone dome to its edge, in the PDB representations. 1036 1037 @keyword mol: The molecule container. 1038 @type mol: MolContainer instance 1039 @keyword cone: The cone object. This should provide the limit_check() method with determines the limits of the distribution accepting two arguments, the polar angle phi and the azimuthal angle theta, and return True if the point is in the limits or False if outside. It should also provide the theta_max() method for returning the theta value for the given phi. 1040 @type cone: class instance 1041 @keyword chain_id: The chain identifier. 1042 @type chain_id: str 1043 @keyword dome_start: The starting atom number of the cone dome residue. 1044 @type dome_start: int 1045 @keyword edge_start: The starting atom number of the cone edge residue. 1046 @type edge_start: int 1047 @keyword scale: The scaling factor to stretch all points by. 1048 @type scale: float 1049 @keyword inc: The number of increments or number of vectors used to generate the outer edge of the cone. 1050 @type inc: int 1051 @keyword distribution: The type of point distribution to use. This can be 'uniform' or 'regular'. 1052 @type distribution: str 1053 """ 1054 1055 # Get the polar and azimuthal angles for the distribution. 1056 if distribution == 'uniform': 1057 phi, theta = angles_uniform(inc) 1058 else: 1059 phi, theta = angles_regular(inc) 1060 1061 # Determine the maximum phi values and the indices of the point just above the edge. 1062 phi_max = zeros(len(theta), float64) 1063 edge_index = zeros(len(theta), int) 1064 for i in range(len(theta)): 1065 # Get the polar angle for the longitude edge atom. 1066 phi_max[i] = cone.phi_max(theta[i]) 1067 1068 # The index. 1069 for j in range(len(phi)): 1070 if phi[j] <= phi_max[i]: 1071 edge_index[i] = j 1072 break 1073 1074 # Reverse edge index. 1075 edge_index_rev = len(phi) - 1 - edge_index 1076 1077 # Debugging. 1078 if debug: 1079 print("\nDome start: %s" % dome_start) 1080 print("Edge start: %s" % edge_start) 1081 print("Edge indices: %s" % edge_index) 1082 print("Edge indices rev: %s" % edge_index_rev) 1083 1084 # Move around the azimuth. 1085 dome_atom = dome_start 1086 edge_atom = edge_start 1087 for i in range(len(theta)): 1088 # Debugging. 1089 if debug: 1090 print("i: %s; theta: %s" % (i, theta[i])) 1091 print("%sDome atom: %s" % (" "*4, get_proton_name(dome_atom))) 1092 print("%sStitching longitudinal line to edge - %s to %s." % (" "*4, get_proton_name(edge_atom), get_proton_name(dome_atom))) 1093 1094 # Connect the two atoms (to stitch up the 2 objects). 1095 mol.atom_connect(index1=dome_atom-1, index2=edge_atom-1) 1096 1097 # Update the edge atom. 1098 edge_atom = edge_atom + 1 1099 1100 # Stitch up the latitude atoms. 1101 for j in range(len(phi)): 1102 # The index for the direction top to bottom! 1103 k = len(phi) - 1 - j 1104 1105 # Debugging. 1106 if debug: 1107 print("%sj: %s; phi: %-20s; k: %s; phi: %-20s; phi_max: %-20s" % (" "*4, j, phi[j], k, phi[k], phi_max[i])) 1108 1109 # No edge. 1110 skip = True 1111 1112 # Forward edge (skip when the latitude is phi max). 1113 fwd_index = i+1 1114 if i == len(theta)-1: 1115 fwd_index = 0 1116 if j >= edge_index[i] and j < edge_index[fwd_index] and not abs(phi_max[fwd_index] - phi[j]) < 1e-6: 1117 # Debugging. 1118 if debug: 1119 print("%sForward edge." % (" "*8)) 1120 1121 # Edge found. 1122 skip = False 1123 forward = True 1124 1125 # Back edge (skip when the latitude is phi max). 1126 rev_index = i-1 1127 if i == 0: 1128 rev_index = len(theta)-1 1129 if i < len(theta)-1 and j > edge_index_rev[i] and j <= edge_index_rev[i+1] and not abs(phi_max[fwd_index] - phi[k]) < 1e-6: 1130 # Debugging. 1131 if debug: 1132 print("%sBack edge." % (" "*8)) 1133 1134 # Edge found. 1135 skip = False 1136 forward = False 1137 1138 # Skipping. 1139 if skip: 1140 continue 1141 1142 # The dome atom to stitch to (current dome atom + one latitude line to shift across). 1143 if forward: 1144 atom = dome_atom + j - edge_index[i] 1145 else: 1146 atom = dome_atom + (edge_index_rev[i]+1) + k - edge_index[fwd_index] 1147 1148 # Debugging. 1149 if debug: 1150 print("%sStitching latitude line to edge - %s to %s." % (" "*8, get_proton_name(edge_atom), get_proton_name(atom))) 1151 1152 # Connect the two atoms (to stitch up the 2 objects). 1153 mol.atom_connect(index1=atom-1, index2=edge_atom-1) 1154 1155 # Increment the cone edge atom number. 1156 edge_atom = edge_atom + 1 1157 1158 # Update the cone dome atom. 1159 dome_atom = dome_atom + (len(phi) - edge_index[i])

1160 1161

1162 -def vect_dist_spherical_angles(inc=20, distribution='uniform'):

1163 """Create a distribution of vectors on a sphere using a distribution of spherical angles. 1164 1165 This function returns an array of unit vectors distributed within 3D space. The unit vectors are generated using the equation:: 1166 1167 | cos(theta) * sin(phi) | 1168 vector = | sin(theta) * sin(phi) |. 1169 | cos(phi) | 1170 1171 The vectors of this distribution generate both longitudinal and latitudinal lines. 1172 1173 1174 @keyword inc: The number of increments in the distribution. 1175 @type inc: int 1176 @keyword distribution: The type of point distribution to use. This can be 'uniform' or 'regular'. 1177 @type distribution: str 1178 @return: The distribution of vectors on a sphere. 1179 @rtype: list of rank-1, 3D numpy arrays 1180 """ 1181 1182 # Get the polar and azimuthal angles for the distribution. 1183 if distribution == 'uniform': 1184 phi, theta = angles_uniform(inc) 1185 else: 1186 phi, theta = angles_regular(inc) 1187 1188 # Initialise array of the distribution of vectors. 1189 vectors = [] 1190 1191 # Loop over the longitudinal lines. 1192 for j in range(len(phi)): 1193 # Loop over the latitudinal lines. 1194 for i in range(len(theta)): 1195 # X coordinate. 1196 x = cos(theta[i]) * sin(phi[j]) 1197 1198 # Y coordinate. 1199 y = sin(theta[i])* sin(phi[j]) 1200 1201 # Z coordinate. 1202 z = cos(phi[j]) 1203 1204 # Append the vector. 1205 vectors.append(array([x, y, z], float64)) 1206 1207 # Return the array of vectors and angles. 1208 return vectors

1209

Source Code for Module generic_fns.structure.geometric