generic_fns.structure.geometric

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2011 Edward d'Auvergne # 4 # # 5 # This file is part of the program relax. # 6 # # 7 # relax 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 2 of the License, or # 10 # (at your option) any later version. # 11 # # 12 # relax 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 relax; if not, write to the Free Software # 19 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA # 20 # # 21 ############################################################################### 22 23 # Python module imports. 24 from math import cos, pi, sin 25 from numpy import arccos, array, dot, eye, float64, transpose, zeros 26 from os import getcwd 27 from string import ascii_uppercase 28 from warnings import warn 29 30 # relax module imports. 31 from generic_fns.mol_res_spin import exists_mol_res_spin_data, spin_loop 32 from generic_fns import pipes 33 from generic_fns.structure.mass import centre_of_mass 34 from internal import Internal 35 from maths_fns.rotation_matrix import two_vect_to_R 36 from relax_errors import RelaxError, RelaxNoPdbError, RelaxNoSequenceError, RelaxNoTensorError, RelaxNoVectorsError 37 from relax_io import get_file_path, open_write_file 38 from relax_warnings import RelaxWarning 39 from status import Status; status = Status() 40 41

42 -def angles_regular(inc=None):

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

73 74

75 -def angles_uniform(inc=None):

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

105 106

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

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

125 126

127 -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):

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

314 315

316 -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):

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

399 400

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

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

582 583

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

585 """Create a PDB representation of the XH vector distribution. 586 587 @param length: The length to set the vectors to in the PDB file. 588 @type length: float 589 @param symmetry: The symmetry flag which if set will create a second PDB chain 'B' which 590 is the same as chain 'A' but with the XH vectors reversed. 591 @type symmetry: bool 592 @param file: The name of the PDB file to create. 593 @type file: str 594 @param dir: The name of the directory to place the PDB file into. 595 @type dir: str 596 @param force: Flag which if set will overwrite any pre-existing file. 597 @type force: bool 598 """ 599 600 # Test if the current pipe exists. 601 pipes.test() 602 603 # Test if a structure has been loaded. 604 if not hasattr(cdp, 'structure') or not cdp.structure.num_models() > 0: 605 raise RelaxNoPdbError 606 607 # Test if sequence data is loaded. 608 if not exists_mol_res_spin_data(): 609 raise RelaxNoSequenceError 610 611 # Test if unit vectors exist. 612 vectors = 0 613 for spin in spin_loop(): 614 if hasattr(spin, 'xh_vect'): 615 vectors = 1 616 break 617 if not vectors: 618 raise RelaxNoVectorsError 619 620 621 # Initialise. 622 ############# 623 624 # Create the structural object. 625 structure = Internal() 626 627 # Add a structure. 628 structure.add_molecule(name='vector_dist') 629 630 # Alias the single molecule from the single model. 631 mol = structure.structural_data[0].mol[0] 632 633 # Initialise the residue and atom numbers. 634 res_num = 1 635 atom_num = 1 636 637 638 # Centre of mass. 639 ################# 640 641 # Calculate the centre of mass. 642 R = centre_of_mass() 643 644 # Increment the residue number. 645 res_num = res_num + 1 646 647 648 # The XH vectors. 649 ################# 650 651 # Loop over the spin systems. 652 for spin, mol_name, res_num, res_name in spin_loop(full_info=True): 653 # Skip deselected spin systems. 654 if not spin.select: 655 continue 656 657 # Skip spin systems missing the xh_vect structure. 658 if not hasattr(spin, 'xh_vect'): 659 continue 660 661 # Scale the vector. 662 vector = spin.xh_vect * length * 1e10 663 664 # Add the central X atom. 665 mol.atom_add(pdb_record='ATOM', atom_num=atom_num, atom_name=spin.name, res_name=res_name, chain_id='A', res_num=res_num, pos=R, segment_id=None, element=spin.element) 666 667 # Add the H atom. 668 mol.atom_add(pdb_record='ATOM', atom_num=atom_num+1, atom_name='H', res_name=res_name, chain_id='A', res_num=res_num, pos=R+vector, segment_id=None, element='H') 669 670 # Connect the two atoms. 671 mol.atom_connect(index1=atom_num-1, index2=atom_num) 672 673 # Increment the atom number. 674 atom_num = atom_num + 2 675 676 # Symmetry chain. 677 if symmetry: 678 # Loop over the spin systems. 679 for spin, mol_name, res_num, res_name in spin_loop(full_info=True): 680 # Skip deselected spin systems. 681 if not spin.select: 682 continue 683 684 # Skip spin systems missing the xh_vect structure. 685 if not hasattr(spin, 'xh_vect'): 686 continue 687 688 # Scale the vector. 689 vector = spin.xh_vect * length * 1e10 690 691 # Add the central X atom. 692 mol.atom_add(pdb_record='ATOM', atom_num=atom_num, atom_name=spin.name, res_name=res_name, chain_id='B', res_num=res_num, pos=R, segment_id=None, element=spin.element) 693 694 # Add the H atom. 695 mol.atom_add(pdb_record='ATOM', atom_num=atom_num+1, atom_name='H', res_name=res_name, chain_id='B', res_num=res_num, pos=R-vector, segment_id=None, element='H') 696 697 # Connect the two atoms. 698 mol.atom_connect(index1=atom_num-1, index2=atom_num) 699 700 # Increment the atom number. 701 atom_num = atom_num + 2 702 703 704 # Create the PDB file. 705 ###################### 706 707 # Print out. 708 print("\nGenerating the PDB file.") 709 710 # Open the PDB file for writing. 711 tensor_pdb_file = open_write_file(file, dir, force=force) 712 713 # Write the data. 714 structure.write_pdb(tensor_pdb_file) 715 716 # Close the file. 717 tensor_pdb_file.close() 718 719 # Add the file to the results file list. 720 if not hasattr(cdp, 'result_files'): 721 cdp.result_files = [] 722 if dir == None: 723 dir = getcwd() 724 cdp.result_files.append(['vector_dist_pdb', 'Vector distribution PDB', get_file_path(file, dir)]) 725 status.observers.result_file.notify()

726 727

728 -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):

729 """Generate a uniformly distributed distribution of atoms on a warped sphere. 730 731 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. 732 733 734 @keyword mol: The molecule container. 735 @type mol: MolContainer instance 736 @keyword res_name: The residue name. 737 @type res_name: str 738 @keyword res_num: The residue number. 739 @type res_num: int 740 @keyword chain_id: The chain identifier. 741 @type chain_id: str 742 @keyword centre: The centre of the distribution. 743 @type centre: numpy array, len 3 744 @keyword R: The optional 3x3 rotation matrix. 745 @type R: 3x3 numpy array 746 @keyword warp: The optional 3x3 warping matrix. 747 @type warp: 3x3 numpy array 748 @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. 749 @type limit_check: function 750 @keyword scale: The scaling factor to stretch all rotated and warped vectors by. 751 @type scale: float 752 @keyword inc: The number of increments or number of vectors used to generate the outer edge of the cone. 753 @type inc: int 754 @keyword distribution: The type of point distribution to use. This can be 'uniform' or 'regular'. 755 @type distribution: str 756 """ 757 758 # Initial atom number. 759 if len(mol.atom_num) == 0: 760 origin_num = 1 761 else: 762 origin_num = mol.atom_num[-1]+1 763 atom_num = origin_num 764 765 # Get the uniform vector distribution. 766 print(" Creating the uniform vector distribution.") 767 vectors = vect_dist_spherical_angles(inc=inc, distribution=distribution) 768 769 # Get the polar and azimuthal angles for the distribution. 770 if distribution == 'uniform': 771 phi, theta = angles_uniform(inc) 772 else: 773 phi, theta = angles_regular(inc) 774 775 # Init the arrays for stitching together. 776 edge = zeros(len(theta)) 777 edge_index = zeros(len(theta), int) 778 edge_phi = zeros(len(theta), float64) 779 edge_atom = zeros(len(theta)) 780 781 # Loop over the radial array of vectors (change in longitude). 782 for i in range(len(theta)): 783 # Debugging. 784 if debug: 785 print("i: %s; theta: %s" % (i, theta[i])) 786 787 # Loop over the vectors of the radial array (change in latitude). 788 for j in range(len(phi)): 789 # Debugging. 790 if debug: 791 print("%sj: %s; phi: %s" % (" "*4, j, phi[j])) 792 793 # Skip the vector if the point is outside of the limits. 794 if limit_check and not limit_check(phi[j], theta[i]): 795 if debug: 796 print("%sOut of cone." % (" "*8)) 797 continue 798 799 # Update the edge for this longitude. 800 if not edge[i]: 801 edge[i] = 1 802 edge_index[i] = j 803 edge_phi[i] = phi[j] 804 edge_atom[i] = atom_num 805 806 # Debugging. 807 if debug: 808 print("%sEdge detected." % (" "*8)) 809 print("%sEdge index: %s" % (" "*8, edge_index[i])) 810 print("%sEdge phi pos: %s" % (" "*8, edge_phi[i])) 811 print("%sEdge atom: %s" % (" "*8, edge_atom[i])) 812 813 # Rotate the vector into the diffusion frame. 814 vector = dot(R, vectors[i + j*len(theta)]) 815 816 # Set the length of the vector to its diffusion rate within the diffusion tensor geometric object. 817 vector = dot(warp, vector) 818 819 # Scale the vector. 820 vector = vector * scale 821 822 # Position relative to the centre of mass. 823 pos = centre + vector 824 825 # Debugging. 826 if debug: 827 print("%sAdding atom %s." % (" "*8, get_proton_name(atom_num))) 828 829 # Add the vector as a H atom of the TNS residue. 830 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') 831 832 # Connect to the previous atom to generate the longitudinal lines (except for the first point). 833 if j > edge_index[i]: 834 # Debugging. 835 if debug: 836 print("%sLongitude line, connecting %s to %s" % (" "*8, get_proton_name(atom_num), get_proton_name(atom_num-1))) 837 838 mol.atom_connect(index1=atom_num-1, index2=atom_num-2) 839 840 # Connect across the radial arrays (to generate the latitudinal lines). 841 if i != 0 and j >= edge_index[i-1]: 842 # The number of atoms back to the previous longitude. 843 num = len(phi) - edge_index[i] 844 845 # Debugging. 846 if debug: 847 print("%sLatitude line, connecting %s to %s" % (" "*8, get_proton_name(atom_num), get_proton_name(atom_num-num))) 848 849 mol.atom_connect(index1=atom_num-1, index2=atom_num-1-num) 850 851 # Connect the last radial array to the first (to zip up the geometric object and close the latitudinal lines). 852 if i == len(theta)-1 and j >= edge_index[0]: 853 # The number of atom in the first longitude line. 854 num = origin_num + j - edge_index[0] 855 856 # Debugging. 857 if debug: 858 print("%sZipping up, connecting %s to %s" % (" "*8, get_proton_name(atom_num), get_proton_name(num))) 859 860 mol.atom_connect(index1=atom_num-1, index2=num-1) 861 862 # Increment the atom number. 863 atom_num = atom_num + 1

864 865

866 -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):

867 """Generate residue representations for the vector and the MC simulationed vectors. 868 869 This is used to create a PDB representation of any vector, including its Monte Carlo 870 simulations. 871 872 @keyword mol: The molecule container. 873 @type mol: MolContainer instance 874 @param vector: The vector to be represented in the PDB. 875 @type vector: numpy array, len 3 876 @param atom_name: The atom name used to label the atom representing the head of the vector. 877 @type atom_name: str 878 @param res_name_vect: The 3 letter PDB residue code used to represent the vector. 879 @type res_name_vect: str 880 @param sim_vectors: The optional Monte Carlo simulation vectors to be represented in the PDB. 881 @type sim_vectors: list of numpy array, each len 3 882 @param res_name_sim: The 3 letter PDB residue code used to represent the Monte Carlo simulation vectors. 883 @type res_name_sim: str 884 @param chain_id: The chain identification code. 885 @type chain_id: str 886 @param res_num: The residue number. 887 @type res_num: int 888 @param origin: The origin for the axis. 889 @type origin: numpy array, len 3 890 @param scale: The scaling factor to stretch the vectors by. 891 @type scale: float 892 @param 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. 893 @type label_placement: float 894 @param neg: If True, then the negative vector positioned at the origin will also be included. 895 @type neg: bool 896 @return: The new residue number. 897 @rtype: int 898 """ 899 900 # The atom numbers (and indices). 901 origin_num = mol.atom_num[-1]+1 902 atom_num = mol.atom_num[-1]+2 903 atom_neg_num = mol.atom_num[-1]+3 904 905 # The origin atom. 906 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') 907 908 # Create the PDB residue representing the vector. 909 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') 910 mol.atom_connect(index1=atom_num-1, index2=origin_num-1) 911 num = 1 912 if neg: 913 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') 914 mol.atom_connect(index1=atom_neg_num-1, index2=origin_num-1) 915 num = 2 916 917 # Add another atom to allow the axis labels to be shifted just outside of the vector itself. 918 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') 919 if neg: 920 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') 921 922 # Print out. 923 print((" " + atom_name + " vector (scaled + shifted to origin): " + repr(origin+vector*scale))) 924 print(" Creating the MC simulation vectors.") 925 926 # Monte Carlo simulations. 927 if sim_vectors != None: 928 for i in range(len(sim_vectors)): 929 # Increment the residue number, so each simulation is a new residue. 930 res_num = res_num + 1 931 932 # The atom numbers (and indices). 933 atom_num = mol.atom_num[-1]+1 934 atom_neg_num = mol.atom_num[-1]+2 935 936 # Create the PDB residue representing the vector. 937 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') 938 mol.atom_connect(index1=atom_num-1, index2=origin_num-1) 939 if neg: 940 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') 941 mol.atom_connect(index1=atom_neg_num-1, index2=origin_num-1) 942 943 # Return the new residue number. 944 return res_num

945 946

947 -def get_proton_name(atom_num):

948 """Return a valid PDB atom name of <4 characters. 949 950 @param atom_num: The number of the atom. 951 @type atom_num: int 952 @return: The atom name to use in the PDB. 953 @rtype: str 954 """ 955 956 # Init the proton first letters and the atom number folding limits. 957 names = ['H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q'] 958 lims = [0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000] 959 960 # Loop over the proton names. 961 for i in range(len(names)): 962 # In the bounds. 963 if atom_num >= lims[i] and atom_num < lims[i+1]: 964 return names[i] + repr(atom_num - lims[i])

965 966

967 -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):

968 """Function for stitching the cone dome to its edge, in the PDB representations. 969 970 @keyword mol: The molecule container. 971 @type mol: MolContainer instance 972 @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. 973 @type cone: class instance 974 @keyword chain_id: The chain identifier. 975 @type chain_id: str 976 @keyword dome_start: The starting atom number of the cone dome residue. 977 @type dome_start: int 978 @keyword edge_start: The starting atom number of the cone edge residue. 979 @type edge_start: int 980 @keyword scale: The scaling factor to stretch all points by. 981 @type scale: float 982 @keyword inc: The number of increments or number of vectors used to generate the outer edge of the cone. 983 @type inc: int 984 @keyword distribution: The type of point distribution to use. This can be 'uniform' or 'regular'. 985 @type distribution: str 986 """ 987 988 # Get the polar and azimuthal angles for the distribution. 989 if distribution == 'uniform': 990 phi, theta = angles_uniform(inc) 991 else: 992 phi, theta = angles_regular(inc) 993 994 # Determine the maximum phi values and the indices of the point just above the edge. 995 phi_max = zeros(len(theta), float64) 996 edge_index = zeros(len(theta), int) 997 for i in range(len(theta)): 998 # Get the polar angle for the longitude edge atom. 999 phi_max[i] = cone.phi_max(theta[i]) 1000 1001 # The index. 1002 for j in range(len(phi)): 1003 if phi[j] <= phi_max[i]: 1004 edge_index[i] = j 1005 break 1006 1007 # Reverse edge index. 1008 edge_index_rev = len(phi) - 1 - edge_index 1009 1010 # Debugging. 1011 if debug: 1012 print("\nDome start: %s" % dome_start) 1013 print("Edge start: %s" % edge_start) 1014 print("Edge indices: %s" % edge_index) 1015 print("Edge indices rev: %s" % edge_index_rev) 1016 1017 # Move around the azimuth. 1018 dome_atom = dome_start 1019 edge_atom = edge_start 1020 for i in range(len(theta)): 1021 # Debugging. 1022 if debug: 1023 print("i: %s; theta: %s" % (i, theta[i])) 1024 print("%sDome atom: %s" % (" "*4, get_proton_name(dome_atom))) 1025 print("%sStitching longitudinal line to edge - %s to %s." % (" "*4, get_proton_name(edge_atom), get_proton_name(dome_atom))) 1026 1027 # Connect the two atoms (to stitch up the 2 objects). 1028 mol.atom_connect(index1=dome_atom-1, index2=edge_atom-1) 1029 1030 # Update the edge atom. 1031 edge_atom = edge_atom + 1 1032 1033 # Stitch up the latitude atoms. 1034 for j in range(len(phi)): 1035 # The index for the direction top to bottom! 1036 k = len(phi) - 1 - j 1037 1038 # Debugging. 1039 if debug: 1040 print("%sj: %s; phi: %-20s; k: %s; phi: %-20s; phi_max: %-20s" % (" "*4, j, phi[j], k, phi[k], phi_max[i])) 1041 1042 # No edge. 1043 skip = True 1044 1045 # Forward edge (skip when the latitude is phi max). 1046 fwd_index = i+1 1047 if i == len(theta)-1: 1048 fwd_index = 0 1049 if j >= edge_index[i] and j < edge_index[fwd_index] and not abs(phi_max[fwd_index] - phi[j]) < 1e-6: 1050 # Debugging. 1051 if debug: 1052 print("%sForward edge." % (" "*8)) 1053 1054 # Edge found. 1055 skip = False 1056 forward = True 1057 1058 # Back edge (skip when the latitude is phi max). 1059 rev_index = i-1 1060 if i == 0: 1061 rev_index = len(theta)-1 1062 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: 1063 # Debugging. 1064 if debug: 1065 print("%sBack edge." % (" "*8)) 1066 1067 # Edge found. 1068 skip = False 1069 forward = False 1070 1071 # Skipping. 1072 if skip: 1073 continue 1074 1075 # The dome atom to stitch to (current dome atom + one latitude line to shift across). 1076 if forward: 1077 atom = dome_atom + j - edge_index[i] 1078 else: 1079 atom = dome_atom + (edge_index_rev[i]+1) + k - edge_index[fwd_index] 1080 1081 # Debugging. 1082 if debug: 1083 print("%sStitching latitude line to edge - %s to %s." % (" "*8, get_proton_name(edge_atom), get_proton_name(atom))) 1084 1085 # Connect the two atoms (to stitch up the 2 objects). 1086 mol.atom_connect(index1=atom-1, index2=edge_atom-1) 1087 1088 # Increment the cone edge atom number. 1089 edge_atom = edge_atom + 1 1090 1091 # Update the cone dome atom. 1092 dome_atom = dome_atom + (len(phi) - edge_index[i])

1093 1094

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

1096 """Create a distribution of vectors on a sphere using a distribution of spherical angles. 1097 1098 This function returns an array of unit vectors distributed within 3D space. The unit vectors are generated using the equation:: 1099 1100 | cos(theta) * sin(phi) | 1101 vector = | sin(theta) * sin(phi) |. 1102 | cos(phi) | 1103 1104 The vectors of this distribution generate both longitudinal and latitudinal lines. 1105 1106 1107 @keyword inc: The number of increments in the distribution. 1108 @type inc: int 1109 @keyword distribution: The type of point distribution to use. This can be 'uniform' or 'regular'. 1110 @type distribution: str 1111 @return: The distribution of vectors on a sphere. 1112 @rtype: list of rank-1, 3D numpy arrays 1113 """ 1114 1115 # Get the polar and azimuthal angles for the distribution. 1116 if distribution == 'uniform': 1117 phi, theta = angles_uniform(inc) 1118 else: 1119 phi, theta = angles_regular(inc) 1120 1121 # Initialise array of the distribution of vectors. 1122 vectors = [] 1123 1124 # Loop over the longitudinal lines. 1125 for j in xrange(len(phi)): 1126 # Loop over the latitudinal lines. 1127 for i in xrange(len(theta)): 1128 # X coordinate. 1129 x = cos(theta[i]) * sin(phi[j]) 1130 1131 # Y coordinate. 1132 y = sin(theta[i])* sin(phi[j]) 1133 1134 # Z coordinate. 1135 z = cos(phi[j]) 1136 1137 # Append the vector. 1138 vectors.append(array([x, y, z], float64)) 1139 1140 # Return the array of vectors and angles. 1141 return vectors

1142

Source Code for Module generic_fns.structure.geometric