pipe_control.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, zeros 25 from os import getcwd 26 27 # relax module imports. 28 from lib.errors import RelaxNoPdbError, RelaxNoSequenceError, RelaxNoVectorsError 29 from lib.io import get_file_path, open_write_file 30 from lib.structure.internal.object import Internal 31 from lib.structure.represent.rotor import rotor_pdb 32 from pipe_control import pipes 33 from pipe_control.interatomic import interatomic_loop 34 from pipe_control.mol_res_spin import exists_mol_res_spin_data, return_spin 35 from pipe_control.structure.mass import pipe_centre_of_mass 36 from status import Status; status = Status() 37 38

39 -def angles_regular(inc=None):

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

70 71

72 -def angles_uniform(inc=None):

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

102 103

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

105 """Create a PDB representation of a rotor motional model. 106 107 @keyword file: The name of the PDB file to create. 108 @type file: str 109 @keyword dir: The name of the directory to place the PDB file into. 110 @type dir: str 111 @keyword rotor_angle: The angle of the rotor motion in degrees. 112 @type rotor_angle: float 113 @keyword axis: The vector defining the rotor axis. 114 @type axis: numpy rank-1, 3D array 115 @keyword axis_pt: A point lying anywhere on the rotor axis. This is used to define the position of the axis in 3D space. 116 @type axis_pt: numpy rank-1, 3D array 117 @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. 118 @type centre: numpy rank-1, 3D array 119 @keyword span: The distance from the central point to the rotor blades (meters). 120 @type span: float 121 @keyword blade_length: The length of the representative rotor blades. 122 @type blade_length: float 123 @keyword force: A flag which if set will overwrite any pre-existing file. 124 @type force: bool 125 @keyword staggered: A flag which if True will cause the rotor blades to be staggered. This is used to avoid blade overlap. 126 @type staggered: bool 127 """ 128 129 # Test if the current pipe exists. 130 pipes.test() 131 132 # Convert the angle to radians. 133 rotor_angle = rotor_angle / 360.0 * 2.0 * pi 134 135 # Create the structural object. 136 structure = Internal() 137 138 # Generate the rotor object. 139 rotor_pdb(structure=structure, rotor_angle=rotor_angle, axis=axis, axis_pt=axis_pt, centre=centre, span=span, blade_length=blade_length, staggered=staggered) 140 141 # Print out. 142 print("\nGenerating the PDB file.") 143 144 # Open the PDB file for writing. 145 tensor_pdb_file = open_write_file(file, dir, force=force) 146 147 # Write the data. 148 structure.write_pdb(tensor_pdb_file) 149 150 # Close the file. 151 tensor_pdb_file.close() 152 153 # Add the file to the results file list. 154 if not hasattr(cdp, 'result_files'): 155 cdp.result_files = [] 156 if dir == None: 157 dir = getcwd() 158 cdp.result_files.append(['rotor_pdb', 'Rotor PDB', get_file_path(file, dir)]) 159 status.observers.result_file.notify()

160 161

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

163 """Create a PDB representation of the vector distribution. 164 165 @keyword length: The length to set the vectors to in the PDB file. 166 @type length: float 167 @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. 168 @type symmetry: bool 169 @keyword file: The name of the PDB file to create. 170 @type file: str 171 @keyword dir: The name of the directory to place the PDB file into. 172 @type dir: str 173 @keyword force: Flag which if set will overwrite any pre-existing file. 174 @type force: bool 175 """ 176 177 # Test if the current pipe exists. 178 pipes.test() 179 180 # Test if a structure has been loaded. 181 if not hasattr(cdp, 'structure') or not cdp.structure.num_models() > 0: 182 raise RelaxNoPdbError 183 184 # Test if sequence data is loaded. 185 if not exists_mol_res_spin_data(): 186 raise RelaxNoSequenceError 187 188 # Test if unit vectors exist. 189 vectors = False 190 for interatom in interatomic_loop(): 191 if hasattr(interatom, 'vector'): 192 vectors = True 193 break 194 if not vectors: 195 raise RelaxNoVectorsError 196 197 198 # Initialise. 199 ############# 200 201 # Create the structural object. 202 structure = Internal() 203 204 # Add a structure. 205 structure.add_molecule(name='vector_dist') 206 207 # Alias the single molecule from the single model. 208 mol = structure.structural_data[0].mol[0] 209 210 # Initialise the residue and atom numbers. 211 res_num = 1 212 atom_num = 1 213 214 215 # Centre of mass. 216 ################# 217 218 # Calculate the centre of mass. 219 R = pipe_centre_of_mass() 220 221 # Increment the residue number. 222 res_num = res_num + 1 223 224 225 # The vectors. 226 ############## 227 228 # Loop over the interatomic data containers. 229 for interatom in interatomic_loop(): 230 # Get the spins. 231 spin1 = return_spin(interatom.spin_id1) 232 spin2 = return_spin(interatom.spin_id2) 233 234 # Skip deselected spin systems. 235 if not spin1.select or not spin2.select: 236 continue 237 238 # Skip containers missing vectors. 239 if not hasattr(interatom, 'vector'): 240 continue 241 242 # Scale the vector. 243 vector = interatom.vector * length * 1e10 244 245 # Add the first spin as the central atom. 246 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) 247 248 # Add the second spin as the end atom. 249 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) 250 251 # Connect the two atoms. 252 mol.atom_connect(index1=atom_num-1, index2=atom_num) 253 254 # Increment the atom number. 255 atom_num = atom_num + 2 256 257 # Symmetry chain. 258 if symmetry: 259 # Loop over the interatomic data containers. 260 for interatom in interatomic_loop(): 261 # Get the spins. 262 spin1 = return_spin(interatom.spin_id1) 263 spin2 = return_spin(interatom.spin_id2) 264 265 # Skip deselected spin systems. 266 if not spin1.select or not spin2.select: 267 continue 268 269 # Skip containers missing vectors. 270 if not hasattr(interatom, 'vector'): 271 continue 272 273 # Scale the vector. 274 vector = interatom.vector * length * 1e10 275 276 # Add the first spin as the central atom. 277 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) 278 279 # Add the second spin as the end atom. 280 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) 281 282 # Connect the two atoms. 283 mol.atom_connect(index1=atom_num-1, index2=atom_num) 284 285 # Increment the atom number. 286 atom_num = atom_num + 2 287 288 289 # Create the PDB file. 290 ###################### 291 292 # Print out. 293 print("\nGenerating the PDB file.") 294 295 # Open the PDB file for writing. 296 tensor_pdb_file = open_write_file(file, dir, force=force) 297 298 # Write the data. 299 structure.write_pdb(tensor_pdb_file) 300 301 # Close the file. 302 tensor_pdb_file.close() 303 304 # Add the file to the results file list. 305 if not hasattr(cdp, 'result_files'): 306 cdp.result_files = [] 307 if dir == None: 308 dir = getcwd() 309 cdp.result_files.append(['vector_dist_pdb', 'Vector distribution PDB', get_file_path(file, dir)]) 310 status.observers.result_file.notify()

311 312

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

314 """Generate a uniformly distributed distribution of atoms on a warped sphere. 315 316 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. 317 318 319 @keyword mol: The molecule container. 320 @type mol: MolContainer instance 321 @keyword res_name: The residue name. 322 @type res_name: str 323 @keyword res_num: The residue number. 324 @type res_num: int 325 @keyword chain_id: The chain identifier. 326 @type chain_id: str 327 @keyword centre: The centre of the distribution. 328 @type centre: numpy array, len 3 329 @keyword R: The optional 3x3 rotation matrix. 330 @type R: 3x3 numpy array 331 @keyword warp: The optional 3x3 warping matrix. 332 @type warp: 3x3 numpy array 333 @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. 334 @type limit_check: function 335 @keyword scale: The scaling factor to stretch all rotated and warped vectors by. 336 @type scale: float 337 @keyword inc: The number of increments or number of vectors. 338 @type inc: int 339 @keyword distribution: The type of point distribution to use. This can be 'uniform' or 'regular'. 340 @type distribution: str 341 """ 342 343 # Initial atom number. 344 if len(mol.atom_num) == 0: 345 origin_num = 1 346 else: 347 origin_num = mol.atom_num[-1]+1 348 atom_num = origin_num 349 350 # Get the uniform vector distribution. 351 print(" Creating the uniform vector distribution.") 352 vectors = vect_dist_spherical_angles(inc=inc, distribution=distribution) 353 354 # Get the polar and azimuthal angles for the distribution. 355 if distribution == 'uniform': 356 phi, theta = angles_uniform(inc) 357 else: 358 phi, theta = angles_regular(inc) 359 360 # Init the arrays for stitching together. 361 edge = zeros(len(theta)) 362 edge_index = zeros(len(theta), int) 363 edge_phi = zeros(len(theta), float64) 364 edge_atom = zeros(len(theta)) 365 366 # Loop over the radial array of vectors (change in longitude). 367 for i in range(len(theta)): 368 # Debugging. 369 if debug: 370 print("i: %s; theta: %s" % (i, theta[i])) 371 372 # Loop over the vectors of the radial array (change in latitude). 373 for j in range(len(phi)): 374 # Debugging. 375 if debug: 376 print("%sj: %s; phi: %s" % (" "*4, j, phi[j])) 377 378 # Skip the vector if the point is outside of the limits. 379 if limit_check and not limit_check(phi[j], theta[i]): 380 if debug: 381 print("%sOut of limits." % (" "*8)) 382 continue 383 384 # Update the edge for this longitude. 385 if not edge[i]: 386 edge[i] = 1 387 edge_index[i] = j 388 edge_phi[i] = phi[j] 389 edge_atom[i] = atom_num 390 391 # Debugging. 392 if debug: 393 print("%sEdge detected." % (" "*8)) 394 print("%sEdge index: %s" % (" "*8, edge_index[i])) 395 print("%sEdge phi pos: %s" % (" "*8, edge_phi[i])) 396 print("%sEdge atom: %s" % (" "*8, edge_atom[i])) 397 398 # Rotate the vector into the frame. 399 vector = dot(R, vectors[i + j*len(theta)]) 400 401 # Warp the vector. 402 vector = dot(warp, vector) 403 404 # Scale the vector. 405 vector = vector * scale 406 407 # Position relative to the centre of mass. 408 pos = centre + vector 409 410 # Debugging. 411 if debug: 412 print("%sAdding atom %s." % (" "*8, get_proton_name(atom_num))) 413 414 # Add the vector as a H atom of the TNS residue. 415 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') 416 417 # Connect to the previous atom to generate the longitudinal lines (except for the first point). 418 if j > edge_index[i]: 419 # Debugging. 420 if debug: 421 print("%sLongitude line, connecting %s to %s" % (" "*8, get_proton_name(atom_num), get_proton_name(atom_num-1))) 422 423 mol.atom_connect(index1=atom_num-1, index2=atom_num-2) 424 425 # Connect across the radial arrays (to generate the latitudinal lines). 426 if i != 0 and j >= edge_index[i-1]: 427 # The number of atoms back to the previous longitude. 428 num = len(phi) - edge_index[i] 429 430 # Debugging. 431 if debug: 432 print("%sLatitude line, connecting %s to %s" % (" "*8, get_proton_name(atom_num), get_proton_name(atom_num-num))) 433 434 mol.atom_connect(index1=atom_num-1, index2=atom_num-1-num) 435 436 # Connect the last radial array to the first (to zip up the geometric object and close the latitudinal lines). 437 if i == len(theta)-1 and j >= edge_index[0]: 438 # The number of atom in the first longitude line. 439 num = origin_num + j - edge_index[0] 440 441 # Debugging. 442 if debug: 443 print("%sZipping up, connecting %s to %s" % (" "*8, get_proton_name(atom_num), get_proton_name(num))) 444 445 mol.atom_connect(index1=atom_num-1, index2=num-1) 446 447 # Increment the atom number. 448 atom_num = atom_num + 1

449 450

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

452 """Generate residue representations for the vector and the MC simulationed vectors. 453 454 This is used to create a PDB representation of any vector, including its Monte Carlo simulations. 455 456 @keyword mol: The molecule container. 457 @type mol: MolContainer instance 458 @keyword vector: The vector to be represented in the PDB. 459 @type vector: numpy array, len 3 460 @keyword atom_name: The atom name used to label the atom representing the head of the vector. 461 @type atom_name: str 462 @keyword res_name_vect: The 3 letter PDB residue code used to represent the vector. 463 @type res_name_vect: str 464 @keyword sim_vectors: The optional Monte Carlo simulation vectors to be represented in the PDB. 465 @type sim_vectors: list of numpy array, each len 3 466 @keyword res_name_sim: The 3 letter PDB residue code used to represent the Monte Carlo simulation vectors. 467 @type res_name_sim: str 468 @keyword chain_id: The chain identification code. 469 @type chain_id: str 470 @keyword res_num: The residue number. 471 @type res_num: int 472 @keyword origin: The origin for the axis. 473 @type origin: numpy array, len 3 474 @keyword scale: The scaling factor to stretch the vectors by. 475 @type scale: float 476 @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. 477 @type label_placement: float 478 @keyword neg: If True, then the negative vector positioned at the origin will also be included. 479 @type neg: bool 480 @return: The new residue number. 481 @rtype: int 482 """ 483 484 # The atom numbers (and indices). 485 origin_num = mol.atom_num[-1]+1 486 atom_num = mol.atom_num[-1]+2 487 atom_neg_num = mol.atom_num[-1]+3 488 489 # The origin atom. 490 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') 491 492 # Create the PDB residue representing the vector. 493 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') 494 mol.atom_connect(index1=atom_num-1, index2=origin_num-1) 495 num = 1 496 if neg: 497 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') 498 mol.atom_connect(index1=atom_neg_num-1, index2=origin_num-1) 499 num = 2 500 501 # Add another atom to allow the axis labels to be shifted just outside of the vector itself. 502 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') 503 if neg: 504 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') 505 506 # Print out. 507 print(" " + atom_name + " vector (scaled + shifted to origin): " + repr(origin+vector*scale)) 508 print(" Creating the MC simulation vectors.") 509 510 # Monte Carlo simulations. 511 if sim_vectors != None: 512 for i in range(len(sim_vectors)): 513 # Increment the residue number, so each simulation is a new residue. 514 res_num = res_num + 1 515 516 # The atom numbers (and indices). 517 atom_num = mol.atom_num[-1]+1 518 atom_neg_num = mol.atom_num[-1]+2 519 520 # Create the PDB residue representing the vector. 521 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') 522 mol.atom_connect(index1=atom_num-1, index2=origin_num-1) 523 if neg: 524 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') 525 mol.atom_connect(index1=atom_neg_num-1, index2=origin_num-1) 526 527 # Return the new residue number. 528 return res_num

529 530

531 -def get_proton_name(atom_num):

532 """Return a valid PDB atom name of <4 characters. 533 534 @param atom_num: The number of the atom. 535 @type atom_num: int 536 @return: The atom name to use in the PDB. 537 @rtype: str 538 """ 539 540 # Init the proton first letters and the atom number folding limits. 541 names = ['H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q'] 542 lims = [0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000] 543 544 # Loop over the proton names. 545 for i in range(len(names)): 546 # In the bounds. 547 if atom_num >= lims[i] and atom_num < lims[i+1]: 548 return names[i] + repr(atom_num - lims[i])

549 550

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

552 """Create a distribution of vectors on a sphere using a distribution of spherical angles. 553 554 This function returns an array of unit vectors distributed within 3D space. The unit vectors are generated using the equation:: 555 556 | cos(theta) * sin(phi) | 557 vector = | sin(theta) * sin(phi) |. 558 | cos(phi) | 559 560 The vectors of this distribution generate both longitudinal and latitudinal lines. 561 562 563 @keyword inc: The number of increments in the distribution. 564 @type inc: int 565 @keyword distribution: The type of point distribution to use. This can be 'uniform' or 'regular'. 566 @type distribution: str 567 @return: The distribution of vectors on a sphere. 568 @rtype: list of rank-1, 3D numpy arrays 569 """ 570 571 # Get the polar and azimuthal angles for the distribution. 572 if distribution == 'uniform': 573 phi, theta = angles_uniform(inc) 574 else: 575 phi, theta = angles_regular(inc) 576 577 # Initialise array of the distribution of vectors. 578 vectors = [] 579 580 # Loop over the longitudinal lines. 581 for j in range(len(phi)): 582 # Loop over the latitudinal lines. 583 for i in range(len(theta)): 584 # X coordinate. 585 x = cos(theta[i]) * sin(phi[j]) 586 587 # Y coordinate. 588 y = sin(theta[i])* sin(phi[j]) 589 590 # Z coordinate. 591 z = cos(phi[j]) 592 593 # Append the vector. 594 vectors.append(array([x, y, z], float64)) 595 596 # Return the array of vectors and angles. 597 return vectors

598

Source Code for Module pipe_control.structure.geometric