generic_fns.structure.main

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2011 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 sqrt 24 from minfx.generic import generic_minimise 25 from numpy import array, dot, float64, ndarray, zeros 26 from numpy.linalg import norm 27 from os import F_OK, access 28 from re import search 29 import sys 30 from warnings import warn 31 32 # relax module imports. 33 from generic_fns import molmol, relax_re 34 from generic_fns.mol_res_spin import create_spin, exists_mol_res_spin_data, generate_spin_id, linear_ave, return_molecule, return_residue, return_spin, spin_loop 35 from generic_fns import pipes 36 from generic_fns.structure.api_base import Displacements 37 from generic_fns.structure.internal import Internal 38 from generic_fns.structure.scientific import Scientific_data 39 from generic_fns.structure.superimpose import fit_to_first, fit_to_mean, Pivot_finder 40 from relax_errors import RelaxError, RelaxFileError, RelaxNoPdbError, RelaxNoSequenceError 41 from relax_io import get_file_path, open_write_file, write_data, write_spin_data 42 from relax_warnings import RelaxWarning, RelaxNoPDBFileWarning, RelaxZeroVectorWarning 43 44

45 -def add_atom(mol_name=None, atom_name=None, res_name=None, res_num=None, pos=[None, None, None], element=None, atom_num=None, chain_id=None, segment_id=None, pdb_record=None):

46 """Add a new atom to the structural data object. 47 48 @keyword mol_name: The name of the molecule. 49 @type mol_name: str 50 @keyword atom_name: The atom name, e.g. 'H1'. 51 @type atom_name: str or None 52 @keyword res_name: The residue name. 53 @type res_name: str or None 54 @keyword res_num: The residue number. 55 @type res_num: int or None 56 @keyword pos: The position vector of coordinates. 57 @type pos: list (length = 3) 58 @keyword element: The element symbol. 59 @type element: str or None 60 @keyword atom_num: The atom number. 61 @type atom_num: int or None 62 @keyword chain_id: The chain identifier. 63 @type chain_id: str or None 64 @keyword segment_id: The segment identifier. 65 @type segment_id: str or None 66 @keyword pdb_record: The optional PDB record name, e.g. 'ATOM' or 'HETATM'. 67 @type pdb_record: str or None 68 """ 69 70 # Test if the current data pipe exists. 71 pipes.test() 72 73 # Place the structural object into the relax data store if needed. 74 if not hasattr(cdp, 'structure'): 75 cdp.structure = Internal() 76 77 # Add the atoms. 78 cdp.structure.add_atom(mol_name=mol_name, atom_name=atom_name, res_name=res_name, res_num=res_num, pos=pos, element=element, atom_num=atom_num, chain_id=chain_id, segment_id=segment_id, pdb_record=pdb_record)

79 80

81 -def connect_atom(index1=None, index2=None):

82 """Connect two atoms. 83 84 @keyword index1: The global index of the first atom. 85 @type index1: str 86 @keyword index2: The global index of the first atom. 87 @type index2: str 88 """ 89 90 # Test if the current data pipe exists. 91 pipes.test() 92 93 # Place the structural object into the relax data store if needed. 94 if not hasattr(cdp, 'structure'): 95 cdp.structure = Internal() 96 97 # Add the atoms. 98 cdp.structure.connect_atom(index1=index1, index2=index2)

99 100

101 -def delete():

102 """Simple function for deleting all structural data.""" 103 104 # Test if the current data pipe exists. 105 pipes.test() 106 107 # Run the object method. 108 cdp.structure.delete() 109 110 # Then remove any spin specific structural info. 111 for spin in spin_loop(): 112 # Delete positional information. 113 if hasattr(spin, 'pos'): 114 del spin.pos 115 116 # Delete bond vectors. 117 if hasattr(spin, 'bond_vect'): 118 del spin.bond_vect 119 if hasattr(spin, 'xh_vect'): 120 del spin.xh_vect

121 122

123 -def displacement(model_from=None, model_to=None, atom_id=None, centroid=None):

124 """Calculate the rotational and translational displacement between two structural models. 125 126 @keyword model_from: The optional model number for the starting position of the displacement. 127 @type model_from: int or None 128 @keyword model_to: The optional model number for the ending position of the displacement. 129 @type model_to: int or None 130 @keyword atom_id: The molecule, residue, and atom identifier string. This matches the spin ID string format. 131 @type atom_id: str or None 132 @keyword centroid: An alternative position of the centroid, used for studying pivoted systems. 133 @type centroid: list of float or numpy rank-1, 3D array 134 """ 135 136 # Test if the current data pipe exists. 137 pipes.test() 138 139 # Convert the model_from and model_to args to lists, is supplied. 140 if model_from != None: 141 model_from = [model_from] 142 if model_to != None: 143 model_to = [model_to] 144 145 # Create a list of all models. 146 models = [] 147 for model in cdp.structure.model_loop(): 148 models.append(model.num) 149 150 # Set model_from or model_to to all models if None. 151 if model_from == None: 152 model_from = models 153 if model_to == None: 154 model_to = models 155 156 # Initialise the data structure. 157 if not hasattr(cdp.structure, 'displacements'): 158 cdp.structure.displacements = Displacements() 159 160 # Loop over the starting models. 161 for i in range(len(model_from)): 162 # Assemble the atomic coordinates. 163 coord_from = [] 164 for pos in cdp.structure.atom_loop(atom_id=atom_id, model_num=model_from[i], pos_flag=True): 165 coord_from.append(pos[0]) 166 167 # Loop over the ending models. 168 for j in range(len(model_to)): 169 # Assemble the atomic coordinates. 170 coord_to = [] 171 for pos in cdp.structure.atom_loop(atom_id=atom_id, model_num=model_to[j], pos_flag=True): 172 coord_to.append(pos[0]) 173 174 # Send to the base container for the calculations. 175 cdp.structure.displacements._calculate(model_from=model_from[i], model_to=model_to[j], coord_from=array(coord_from), coord_to=array(coord_to), centroid=centroid)

176 177

178 -def find_pivot(models=None, atom_id=None, init_pos=None):

179 """Superimpose a set of structural models. 180 181 @keyword models: The list of models to use. If set to None, then all models will be used. 182 @type models: list of int or None 183 @keyword atom_id: The molecule, residue, and atom identifier string. This matches the spin ID string format. 184 @type atom_id: str or None 185 @keyword init_pos: The starting pivot position for the pivot point optimisation. 186 @type init_pos: list of float or numpy rank-1, 3D array 187 """ 188 189 # Test if the current data pipe exists. 190 pipes.test() 191 192 # Initialised the starting position if needed. 193 if init_pos == None: 194 init_pos = zeros(3, float64) 195 init_pos = array(init_pos) 196 197 # Validate the models. 198 cdp.structure.validate_models() 199 200 # Create a list of all models. 201 if models == None: 202 models = [] 203 for model in cdp.structure.model_loop(): 204 models.append(model.num) 205 206 # Assemble the atomic coordinates of all models. 207 coord = [] 208 for model in models: 209 coord.append([]) 210 for pos in cdp.structure.atom_loop(atom_id=atom_id, model_num=model, pos_flag=True): 211 coord[-1].append(pos[0]) 212 coord[-1] = array(coord[-1]) 213 coord = array(coord) 214 215 # The target function. 216 finder = Pivot_finder(models, coord) 217 results = generic_minimise(func=finder.func, x0=init_pos, min_algor='simplex', print_flag=1) 218 219 # No result. 220 if results == None: 221 return 222 223 # Store the data. 224 cdp.structure.pivot = results 225 226 # Print out. 227 print("Motional pivot found at: %s" % results)

228 229

230 -def get_pos(spin_id=None, str_id=None, ave_pos=False):

231 """Load the spins from the structural object into the relax data store. 232 233 @keyword spin_id: The molecule, residue, and spin identifier string. 234 @type spin_id: str 235 @keyword str_id: The structure identifier. This can be the file name, model number, or structure number. 236 @type str_id: int or str 237 @keyword ave_pos: A flag specifying if the average atom position or the atom position from all loaded structures is loaded into the SpinContainer. 238 @type ave_pos: bool 239 """ 240 241 # Test if the current data pipe exists. 242 pipes.test() 243 244 # Test if the structure exists. 245 if not hasattr(cdp, 'structure') or not cdp.structure.num_models() or not cdp.structure.num_molecules(): 246 raise RelaxNoPdbError 247 248 # Loop over all atoms of the spin_id selection. 249 model_index = -1 250 last_model = None 251 data = [] 252 for model_num, mol_name, res_num, res_name, atom_num, atom_name, element, pos in cdp.structure.atom_loop(atom_id=spin_id, str_id=str_id, model_num_flag=True, mol_name_flag=True, res_num_flag=True, res_name_flag=True, atom_num_flag=True, atom_name_flag=True, element_flag=True, pos_flag=True, ave=ave_pos): 253 # Update the model info. 254 if last_model != model_num: 255 model_index = model_index + 1 256 last_model = model_num 257 258 # Remove the '+' regular expression character from the mol, res, and spin names! 259 if mol_name and search('\+', mol_name): 260 mol_name = mol_name.replace('+', '') 261 if res_name and search('\+', res_name): 262 res_name = res_name.replace('+', '') 263 if atom_name and search('\+', atom_name): 264 atom_name = atom_name.replace('+', '') 265 266 # The spin identification string. The residue name and spin num is not included to allow molecules with point mutations to be used as different models. 267 id = generate_spin_id(res_num=res_num, res_name=None, spin_name=atom_name) 268 269 # Get the spin container. 270 spin_cont = return_spin(id) 271 272 # Skip the spin if it doesn't exist. 273 if spin_cont == None: 274 continue 275 276 # Add the position vector to the spin container. 277 if ave_pos: 278 spin_cont.pos = pos 279 else: 280 if not hasattr(spin_cont, 'pos'): 281 spin_cont.pos = [] 282 spin_cont.pos.append(pos) 283 284 # Store the data for a printout at the end. 285 data.append([id, repr(pos)]) 286 287 # No positions found. 288 if not len(data): 289 raise RelaxError("No positional information matching the spin ID '%s' could be found." % spin_id) 290 291 # Update pseudo-atoms. 292 for spin in spin_loop(): 293 if hasattr(spin, 'members'): 294 # Get the spin positions. 295 positions = [] 296 for atom in spin.members: 297 # Get the spin container. 298 subspin = return_spin(atom) 299 300 # Test that the spin exists. 301 if subspin == None: 302 raise RelaxNoSpinError(atom) 303 304 # Test the position. 305 if not hasattr(subspin, 'pos') or not subspin.pos: 306 raise RelaxError("Positional information is not available for the atom '%s'." % atom) 307 308 # Alias the position. 309 pos = subspin.pos 310 311 # Convert to a list of lists if not already. 312 multi_model = True 313 if type(pos[0]) in [float, float64]: 314 multi_model = False 315 pos = [pos] 316 317 # Store the position. 318 positions.append([]) 319 for i in range(len(pos)): 320 positions[-1].append(pos[i].tolist()) 321 322 # The averaging. 323 if spin.averaging == 'linear': 324 # Average pos. 325 ave = linear_ave(positions) 326 327 # Convert to the correct structure. 328 if multi_model: 329 spin.pos = ave 330 else: 331 spin.pos = ave[0] 332 333 # Print out. 334 write_data(out=sys.stdout, headings=["Spin_ID", "Position"], data=data)

335 336

337 -def load_spins(spin_id=None, str_id=None, ave_pos=False):

338 """Load the spins from the structural object into the relax data store. 339 340 @keyword spin_id: The molecule, residue, and spin identifier string. 341 @type spin_id: str 342 @keyword str_id: The structure identifier. This can be the file name, model number, or structure number. 343 @type str_id: int or str 344 @keyword ave_pos: A flag specifying if the average atom position or the atom position from all loaded structures is loaded into the SpinContainer. 345 @type ave_pos: bool 346 """ 347 348 # Test if the current data pipe exists. 349 pipes.test() 350 351 # Test if the structure exists. 352 if not hasattr(cdp, 'structure') or not cdp.structure.num_models() or not cdp.structure.num_molecules(): 353 raise RelaxNoPdbError 354 355 # Print out. 356 print("Adding the following spins to the relax data store.\n") 357 358 # Init the data for printing out. 359 mol_names = [] 360 res_nums = [] 361 res_names = [] 362 spin_nums = [] 363 spin_names = [] 364 365 # Initialise data for the atom loop. 366 model_index = -1 367 last_model = 1000000 368 369 # Loop over all atoms of the spin_id selection. 370 for model_num, mol_name, res_num, res_name, atom_num, atom_name, element, pos in cdp.structure.atom_loop(atom_id=spin_id, str_id=str_id, model_num_flag=True, mol_name_flag=True, res_num_flag=True, res_name_flag=True, atom_num_flag=True, atom_name_flag=True, element_flag=True, pos_flag=True, ave=ave_pos): 371 # Update the model info. 372 if last_model != model_num: 373 model_index = model_index + 1 374 last_model = model_num 375 376 # Remove the '+' regular expression character from the mol, res, and spin names! 377 if mol_name and search('\+', mol_name): 378 mol_name = mol_name.replace('+', '') 379 if res_name and search('\+', res_name): 380 res_name = res_name.replace('+', '') 381 if atom_name and search('\+', atom_name): 382 atom_name = atom_name.replace('+', '') 383 384 # Generate a spin ID for the current atom. 385 id = generate_spin_id(mol_name=mol_name, res_num=res_num, res_name=res_name, spin_num=atom_num, spin_name=atom_name) 386 387 # Create the spin. 388 try: 389 spin_cont = create_spin(mol_name=mol_name, res_num=res_num, res_name=res_name, spin_num=atom_num, spin_name=atom_name) 390 except RelaxError: 391 # Throw a warning if still on the first model. 392 if not model_index: 393 warn(RelaxWarning("The spin '%s' already exists." % id)) 394 continue 395 396 # Otherwise, get the spin container. 397 spin_cont = return_spin(id) 398 399 # Append all the spin ID info for the first model for printing later. 400 if model_index == 0: 401 mol_names.append(mol_name) 402 res_nums.append(res_num) 403 res_names.append(res_name) 404 spin_nums.append(atom_num) 405 spin_names.append(atom_name) 406 407 # Convert the position vector to a numpy array. 408 pos = array(pos, float64) 409 410 # Average position vector (already averaged across models in the atom_loop). 411 if ave_pos: 412 spin_cont.pos = pos 413 414 # All positions. 415 else: 416 # Initialise. 417 if not hasattr(spin_cont, 'pos'): 418 spin_cont.pos = [] 419 420 # Add the current model's position. 421 spin_cont.pos.append(pos) 422 423 # Add the element. 424 if not model_index: 425 spin_cont.element = element 426 427 # Catch no data. 428 if len(mol_names) == 0: 429 warn(RelaxWarning("No spins matching the '%s' ID string could be found." % spin_id)) 430 return 431 432 # Print out. 433 write_spin_data(file=sys.stdout, mol_names=mol_names, res_nums=res_nums, res_names=res_names, spin_nums=spin_nums, spin_names=spin_names)

434 435

436 -def read_pdb(file=None, dir=None, read_mol=None, set_mol_name=None, read_model=None, set_model_num=None, parser='internal', verbosity=1, fail=True):

437 """The PDB loading function. 438 439 Parsers 440 ======= 441 442 A number of parsers are available for reading PDB files. These include: 443 444 - 'scientific', the Scientific Python PDB parser. 445 - 'internal', a low quality yet fast PDB parser built into relax. 446 447 448 @keyword file: The name of the PDB file to read. 449 @type file: str 450 @keyword dir: The directory where the PDB file is located. If set to None, then the file will be searched for in the current directory. 451 @type dir: str or None 452 @keyword read_mol: The molecule(s) to read from the file, independent of model. The molecules are determined differently by the different parsers, but are numbered consecutively from 1. If set to None, then all molecules will be loaded. 453 @type read_mol: None, int, or list of int 454 @keyword set_mol_name: Set the names of the molecules which are loaded. If set to None, then the molecules will be automatically labelled based on the file name or other information. 455 @type set_mol_name: None, str, or list of str 456 @keyword read_model: The PDB model to extract from the file. If set to None, then all models will be loaded. 457 @type read_model: None, int, or list of int 458 @keyword set_model_num: Set the model number of the loaded molecule. If set to None, then the PDB model numbers will be preserved, if they exist. 459 @type set_model_num: None, int, or list of int 460 @keyword parser: The parser to be used to read the PDB file. 461 @type parser: str 462 @keyword fail: A flag which, if True, will cause a RelaxError to be raised if the PDB file does not exist. If False, then a RelaxWarning will be trown instead. 463 @type fail: bool 464 @keyword verbosity: The amount of information to print to screen. Zero corresponds to minimal output while higher values increase the amount of output. The default value is 1. 465 @type verbosity: int 466 @raise RelaxFileError: If the fail flag is set, then a RelaxError is raised if the PDB file does not exist. 467 """ 468 469 # Test if the current data pipe exists. 470 pipes.test() 471 472 # The file path. 473 file_path = get_file_path(file, dir) 474 475 # Try adding file extensions to the end of the file path, if the file can't be found. 476 file_path_orig = file_path 477 if not access(file_path, F_OK): 478 # List of possible extensions. 479 for ext in ['.pdb', '.gz', '.pdb.gz', '.bz2', '.pdb.bz2']: 480 # Add the extension if the file can be found. 481 if access(file_path+ext, F_OK): 482 file_path = file_path + ext 483 484 # Test if the file exists. 485 if not access(file_path, F_OK): 486 if fail: 487 raise RelaxFileError('PDB', file_path_orig) 488 else: 489 warn(RelaxNoPDBFileWarning(file_path)) 490 return 491 492 # Check that the parser is the same as the currently loaded PDB files. 493 if hasattr(cdp, 'structure') and cdp.structure.id != parser: 494 raise RelaxError("The " + repr(parser) + " parser does not match the " + repr(cdp.structure.id) + " parser of the PDB loaded into the current pipe.") 495 496 # Place the parser specific structural object into the relax data store. 497 if not hasattr(cdp, 'structure'): 498 if parser == 'scientific': 499 cdp.structure = Scientific_data() 500 elif parser == 'internal': 501 cdp.structure = Internal() 502 503 # Load the structures. 504 cdp.structure.load_pdb(file_path, read_mol=read_mol, set_mol_name=set_mol_name, read_model=read_model, set_model_num=set_model_num, verbosity=verbosity) 505 506 # Load into Molmol (if running). 507 molmol.molmol_obj.open_pdb()

508 509

510 -def read_xyz(file=None, dir=None, read_mol=None, set_mol_name=None, read_model=None, set_model_num=None, verbosity=1, fail=True):

511 """The XYZ loading function. 512 513 514 @keyword file: The name of the XYZ file to read. 515 @type file: str 516 @keyword dir: The directory where the XYZ file is located. If set to None, then the 517 file will be searched for in the current directory. 518 @type dir: str or None 519 @keyword read_mol: The molecule(s) to read from the file, independent of model. 520 If set to None, then all molecules will be loaded. 521 @type read_mol: None, int, or list of int 522 @keyword set_mol_name: Set the names of the molecules which are loaded. If set to None, then 523 the molecules will be automatically labelled based on the file name or 524 other information. 525 @type set_mol_name: None, str, or list of str 526 @keyword read_model: The XYZ model to extract from the file. If set to None, then all models 527 will be loaded. 528 @type read_model: None, int, or list of int 529 @keyword set_model_num: Set the model number of the loaded molecule. If set to None, then the 530 XYZ model numbers will be preserved, if they exist. 531 @type set_model_num: None, int, or list of int 532 @keyword fail: A flag which, if True, will cause a RelaxError to be raised if the XYZ 533 file does not exist. If False, then a RelaxWarning will be trown 534 instead. 535 @type fail: bool 536 @keyword verbosity: The amount of information to print to screen. Zero corresponds to 537 minimal output while higher values increase the amount of output. The 538 default value is 1. 539 @type verbosity: int 540 @raise RelaxFileError: If the fail flag is set, then a RelaxError is raised if the XYZ file 541 does not exist. 542 """ 543 544 # Test if the current data pipe exists. 545 pipes.test() 546 547 # The file path. 548 file_path = get_file_path(file, dir) 549 550 # Try adding '.xyz' to the end of the file path, if the file can't be found. 551 if not access(file_path, F_OK): 552 file_path_orig = file_path 553 file_path = file_path + '.xyz' 554 555 # Test if the file exists. 556 if not access(file_path, F_OK): 557 if fail: 558 raise RelaxFileError('XYZ', file_path_orig) 559 else: 560 warn(RelaxNoPDBFileWarning(file_path)) 561 return 562 563 # Place the structural object into the relax data store. 564 if not hasattr(cdp, 'structure'): 565 cdp.structure = Internal() 566 567 # Load the structures. 568 cdp.structure.load_xyz(file_path, read_mol=read_mol, set_mol_name=set_mol_name, read_model=read_model, set_model_num=set_model_num, verbosity=verbosity)

569 570

571 -def rotate(R=None, origin=None, model=None, atom_id=None):

572 """Rotate the structural data about the origin by the specified forwards rotation. 573 574 @keyword R: The forwards rotation matrix. 575 @type R: numpy 3D, rank-2 array or a 3x3 list of floats 576 @keyword origin: The origin of the rotation. If not supplied, the origin will be set to [0, 0, 0]. 577 @type origin: numpy 3D, rank-1 array or list of len 3 or None 578 @keyword model: The model to rotate. If None, all models will be rotated. 579 @type model: int 580 @keyword atom_id: The molecule, residue, and atom identifier string. Only atoms matching this selection will be used. 581 @type atom_id: str or None 582 """ 583 584 # Test if the current data pipe exists. 585 pipes.test() 586 587 # Test if the structure exists. 588 if not hasattr(cdp, 'structure') or not cdp.structure.num_models() or not cdp.structure.num_molecules(): 589 raise RelaxNoPdbError 590 591 # Set the origin if not supplied. 592 if origin == None: 593 origin = [0, 0, 0] 594 595 # Convert the args to numpy float data structures. 596 R = array(R, float64) 597 origin = array(origin, float64) 598 599 # Call the specific code. 600 cdp.structure.rotate(R=R, origin=origin, model=model, atom_id=atom_id)

601 602

603 -def set_vector(spin=None, xh_vect=None):

604 """Place the XH unit vector into the spin container object. 605 606 @keyword spin: The spin container object. 607 @type spin: SpinContainer instance 608 @keyword xh_vect: The unit vector parallel to the XH bond. 609 @type xh_vect: array of len 3 610 """ 611 612 # Place the XH unit vector into the container. 613 spin.xh_vect = xh_vect

614 615

616 -def superimpose(models=None, method='fit to mean', atom_id=None, centroid=None):

617 """Superimpose a set of structural models. 618 619 @keyword models: The list of models to superimpose. If set to None, then all models will be used. 620 @type models: list of int or None 621 @keyword method: The superimposition method. It must be one of 'fit to mean' or 'fit to first'. 622 @type method: str 623 @keyword atom_id: The molecule, residue, and atom identifier string. This matches the spin ID string format. 624 @type atom_id: str or None 625 @keyword centroid: An alternative position of the centroid to allow for different superpositions, for example of pivot point motions. 626 @type centroid: list of float or numpy rank-1, 3D array 627 """ 628 629 # Check the method. 630 allowed = ['fit to mean', 'fit to first'] 631 if method not in allowed: 632 raise RelaxError("The superimposition method '%s' is unknown. It must be one of %s." % (method, allowed)) 633 634 # Test if the current data pipe exists. 635 pipes.test() 636 637 # Validate the models. 638 cdp.structure.validate_models() 639 640 # Create a list of all models. 641 if models == None: 642 models = [] 643 for model in cdp.structure.model_loop(): 644 models.append(model.num) 645 646 # Assemble the atomic coordinates of all models. 647 coord = [] 648 for model in models: 649 coord.append([]) 650 for pos in cdp.structure.atom_loop(atom_id=atom_id, model_num=model, pos_flag=True): 651 coord[-1].append(pos[0]) 652 coord[-1] = array(coord[-1]) 653 654 # The different algorithms. 655 if method == 'fit to mean': 656 T, R, pivot = fit_to_mean(models=models, coord=coord, centroid=centroid) 657 elif method == 'fit to first': 658 T, R, pivot = fit_to_first(models=models, coord=coord, centroid=centroid) 659 660 661 # Update to the new coordinates. 662 for i in range(len(models)): 663 # Translate the molecule first (the rotational pivot is defined in the first model). 664 translate(T=T[i], model=models[i]) 665 666 # Rotate the molecule. 667 rotate(R=R[i], origin=pivot[i], model=models[i])

668 669

670 -def translate(T=None, model=None, atom_id=None):

671 """Shift the structural data by the specified translation vector. 672 673 @keyword T: The translation vector. 674 @type T: numpy rank-1, 3D array or list of float 675 @keyword model: The model to translate. If None, all models will be rotated. 676 @type model: int or None 677 @keyword atom_id: The molecule, residue, and atom identifier string. Only atoms matching this selection will be used. 678 @type atom_id: str or None 679 """ 680 681 # Test if the current data pipe exists. 682 pipes.test() 683 684 # Test if the structure exists. 685 if not hasattr(cdp, 'structure') or not cdp.structure.num_models() or not cdp.structure.num_molecules(): 686 raise RelaxNoPdbError 687 688 # Convert the args to numpy float data structures. 689 T = array(T, float64) 690 691 # Call the specific code. 692 cdp.structure.translate(T=T, model=model, atom_id=atom_id)

693 694

695 -def vectors(spin_id1=None, spin_id2=None, model=None, verbosity=1, ave=True, unit=True):

696 """Extract the bond vectors from the loaded structures and store them in the spin container. 697 698 @keyword spin_id1: The spin identifier string of the first spin of the pair. 699 @type spin_id1: str 700 @keyword spin_id2: The spin identifier string of the second spin of the pair. 701 @type spin_id2: str 702 @keyword model: The model to extract the vector from. If None, all vectors will be extracted. 703 @type model: str 704 @keyword verbosity: The higher the value, the more information is printed to screen. 705 @type verbosity: int 706 @keyword ave: A flag which if True will cause the average of all vectors to be extracted. 707 @type ave: bool 708 @keyword unit: A flag which if True will cause the function to calculate the unit vectors. 709 @type unit: bool 710 """ 711 712 # Test if the current data pipe exists. 713 pipes.test() 714 715 # Test if the PDB file has been loaded. 716 if not hasattr(cdp, 'structure'): 717 raise RelaxNoPdbError 718 719 # Test if sequence data is loaded. 720 if not exists_mol_res_spin_data(): 721 raise RelaxNoSequenceError 722 723 # Print out. 724 if verbosity: 725 # Number of models. 726 num_models = cdp.structure.num_models() 727 728 # Multiple models loaded. 729 if num_models > 1: 730 if model: 731 print("Extracting vectors for model '%s'." % model) 732 else: 733 print("Extracting vectors for all %s models." % num_models) 734 if ave: 735 print("Averaging all vectors.") 736 737 # Single model loaded. 738 else: 739 print("Extracting vectors from the single model.") 740 741 # Unit vectors. 742 if unit: 743 print("Calculating the unit vectors.") 744 745 # Loop over the spins. 746 no_vectors = True 747 for spin, mol_name, res_num, res_name in spin_loop(selection=spin_id, full_info=True): 748 # Skip deselected spins. 749 if not spin.select: 750 continue 751 752 # The spin identification string. The residue name and spin num is not included to allow molecules with point mutations to be used as different models. 753 id = generate_spin_id(res_num=res_num, res_name=None, spin_name=spin.name, spin_num=spin.num) 754 755 # Test that the spin number or name are set (one or both are essential for the identification of the atom). 756 if spin.num == None and spin.name == None: 757 warn(RelaxWarning("Either the spin number or name must be set for the spin " + repr(id) + " to identify the corresponding atom in the molecule.")) 758 continue 759 760 # The bond vector already exists. 761 if hasattr(spin, 'vector'): 762 obj = getattr(spin, 'vector') 763 if obj != None: 764 warn(RelaxWarning("The bond vector for the spin " + repr(id) + " already exists.")) 765 continue 766 767 # Get the bond info. 768 bond_vectors, attached_name, warnings = cdp.structure.bond_vectors(attached_atom=attached, model_num=model, res_num=res_num, spin_name=spin.name, spin_num=spin.num, return_name=True, return_warnings=True) 769 770 # No attached atom. 771 if not bond_vectors: 772 # Warning messages. 773 if warnings: 774 warn(RelaxWarning(warnings + " (atom ID " + repr(id) + ").")) 775 776 # Skip the spin. 777 continue 778 779 # Set the attached atom name. 780 if not hasattr(spin, 'attached_atom'): 781 spin.attached_atom = attached_name 782 elif spin.attached_atom != attached_name: 783 raise RelaxError("The " + repr(spin.attached_atom) + " atom already attached to the spin does not match the attached atom " + repr(attached_name) + ".") 784 785 # Initialise the average vector. 786 if ave: 787 ave_vector = zeros(3, float64) 788 789 # Loop over the individual vectors. 790 for i in range(len(bond_vectors)): 791 # Unit vector. 792 if unit: 793 # Normalisation factor. 794 norm_factor = norm(bond_vectors[i]) 795 796 # Test for zero length. 797 if norm_factor == 0.0: 798 warn(RelaxZeroVectorWarning(id)) 799 800 # Calculate the normalised vector. 801 else: 802 bond_vectors[i] = bond_vectors[i] / norm_factor 803 804 # Sum the vectors. 805 if ave: 806 ave_vector = ave_vector + bond_vectors[i] 807 808 # Average. 809 if ave: 810 vector = ave_vector / float(len(bond_vectors)) 811 else: 812 vector = bond_vectors 813 814 # Convert to a single vector if needed. 815 if len(vector) == 1: 816 vector = vector[0] 817 818 # Set the vector. 819 setattr(spin, 'vector', vector) 820 821 # We have a vector! 822 no_vectors = False 823 824 # Print out of modified spins. 825 if verbosity: 826 # The number of vectors. 827 num = len(bond_vectors) 828 plural = 's' 829 if num == 1: 830 plural = '' 831 832 if spin.name: 833 print("Extracted %s %s-%s vector%s for the spin '%s'." % (num, spin.name, attached_name, plural, id)) 834 else: 835 print("Extracted %s %s-%s vector%s for the spin '%s'." % (num, spin.num, attached_name, plural, id)) 836 837 # Right, catch the problem of missing vectors to prevent massive user confusion! 838 if no_vectors: 839 raise RelaxError("No vectors could be extracted.")

840 841

842 -def write_pdb(file=None, dir=None, model_num=None, compress_type=0, force=False):

843 """The PDB writing function. 844 845 @keyword file: The name of the PDB file to write. 846 @type file: str 847 @keyword dir: The directory where the PDB file will be placed. If set to None, then the file will be placed in the current directory. 848 @type dir: str or None 849 @keyword model_num: The model to place into the PDB file. If not supplied, then all models will be placed into the file. 850 @type model_num: None or int 851 @keyword compress_type: The compression type. The integer values correspond to the compression type: 0, no compression; 1, Bzip2 compression; 2, Gzip compression. 852 @type compress_type: int 853 @keyword force: The force flag which if True will cause the file to be overwritten. 854 @type force: bool 855 """ 856 857 # Test if the current data pipe exists. 858 pipes.test() 859 860 # Check if the structural object exists. 861 if not hasattr(cdp, 'structure'): 862 raise RelaxError("No structural data is present in the current data pipe.") 863 864 # Check if the structural object is writable. 865 if cdp.structure.id in ['scientific']: 866 raise RelaxError("The structures from the " + cdp.structure.id + " parser are not writable.") 867 868 # The file path. 869 file_path = get_file_path(file, dir) 870 871 # Add '.pdb' to the end of the file path if it isn't there yet. 872 if not search(".pdb$", file_path): 873 file_path = file_path + '.pdb' 874 875 # Open the file for writing. 876 file = open_write_file(file_path, compress_type=compress_type, force=force) 877 878 # Write the structures. 879 cdp.structure.write_pdb(file, model_num=model_num)

880

Source Code for Module generic_fns.structure.main