generic_fns.structure.geometric

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

733 734

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

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

871 872

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

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

952 953

954 -def get_proton_name(atom_num):

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

972 973

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

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

1100 1101

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

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

1149

Source Code for Module generic_fns.structure.geometric