Package generic_fns :: Package structure :: Module main
[hide private]
[frames] | no frames]

Source Code for Module generic_fns.structure.main

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2003-2011 Edward d'Auvergne                                   # 
  4  #                                                                             # 
  5  # This file is part of the program relax.                                     # 
  6  #                                                                             # 
  7  # relax is free software; you can redistribute it and/or modify               # 
  8  # it under the terms of the GNU General Public License as published by        # 
  9  # the Free Software Foundation; either version 2 of the License, or           # 
 10  # (at your option) any later version.                                         # 
 11  #                                                                             # 
 12  # relax is distributed in the hope that it will be useful,                    # 
 13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 15  # GNU General Public License for more details.                                # 
 16  #                                                                             # 
 17  # You should have received a copy of the GNU General Public License           # 
 18  # along with relax; if not, write to the Free Software                        # 
 19  # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA   # 
 20  #                                                                             # 
 21  ############################################################################### 
 22   
 23  # Python module imports. 
 24  from math import sqrt 
 25  from minfx.generic import generic_minimise 
 26  from numpy import array, dot, float64, ndarray, zeros 
 27  from numpy.linalg import norm 
 28  from os import F_OK, access 
 29  from re import search 
 30  from string import replace 
 31  import sys 
 32  from warnings import warn 
 33   
 34  # relax module imports. 
 35  from generic_fns import molmol, relax_re 
 36  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 
 37  from generic_fns import pipes 
 38  from generic_fns.structure.api_base import Displacements 
 39  from generic_fns.structure.internal import Internal 
 40  from generic_fns.structure.scientific import Scientific_data 
 41  from generic_fns.structure.superimpose import fit_to_first, fit_to_mean, Pivot_finder 
 42  from relax_errors import RelaxError, RelaxFileError, RelaxNoPdbError, RelaxNoSequenceError 
 43  from relax_io import get_file_path, open_write_file, write_spin_data 
 44  from relax_warnings import RelaxWarning, RelaxNoPDBFileWarning, RelaxZeroVectorWarning 
 45   
 46   
47 -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):
48 """Add a new atom to the structural data object. 49 50 @keyword mol_name: The name of the molecule. 51 @type mol_name: str 52 @keyword atom_name: The atom name, e.g. 'H1'. 53 @type atom_name: str or None 54 @keyword res_name: The residue name. 55 @type res_name: str or None 56 @keyword res_num: The residue number. 57 @type res_num: int or None 58 @keyword pos: The position vector of coordinates. 59 @type pos: list (length = 3) 60 @keyword element: The element symbol. 61 @type element: str or None 62 @keyword atom_num: The atom number. 63 @type atom_num: int or None 64 @keyword chain_id: The chain identifier. 65 @type chain_id: str or None 66 @keyword segment_id: The segment identifier. 67 @type segment_id: str or None 68 @keyword pdb_record: The optional PDB record name, e.g. 'ATOM' or 'HETATM'. 69 @type pdb_record: str or None 70 """ 71 72 # Test if the current data pipe exists. 73 pipes.test() 74 75 # Place the structural object into the relax data store if needed. 76 if not hasattr(cdp, 'structure'): 77 cdp.structure = Internal() 78 79 # Add the atoms. 80 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)
81 82
83 -def connect_atom(index1=None, index2=None):
84 """Connect two atoms. 85 86 @keyword index1: The global index of the first atom. 87 @type index1: str 88 @keyword index2: The global index of the first atom. 89 @type index2: str 90 """ 91 92 # Test if the current data pipe exists. 93 pipes.test() 94 95 # Place the structural object into the relax data store if needed. 96 if not hasattr(cdp, 'structure'): 97 cdp.structure = Internal() 98 99 # Add the atoms. 100 cdp.structure.connect_atom(index1=index1, index2=index2)
101 102
103 -def delete():
104 """Simple function for deleting all structural data.""" 105 106 # Test if the current data pipe exists. 107 pipes.test() 108 109 # Run the object method. 110 cdp.structure.delete() 111 112 # Then remove any spin specific structural info. 113 for spin in spin_loop(): 114 # Delete positional information. 115 if hasattr(spin, 'pos'): 116 del spin.pos 117 118 # Delete bond vectors. 119 if hasattr(spin, 'bond_vect'): 120 del spin.bond_vect 121 if hasattr(spin, 'xh_vect'): 122 del spin.xh_vect
123 124
125 -def displacement(model_from=None, model_to=None, atom_id=None, centroid=None):
126 """Calculate the rotational and translational displacement between two structural models. 127 128 @keyword model_from: The optional model number for the starting position of the displacement. 129 @type model_from: int or None 130 @keyword model_to: The optional model number for the ending position of the displacement. 131 @type model_to: int or None 132 @keyword atom_id: The molecule, residue, and atom identifier string. This matches the spin ID string format. 133 @type atom_id: str or None 134 @keyword centroid: An alternative position of the centroid, used for studying pivoted systems. 135 @type centroid: list of float or numpy rank-1, 3D array 136 """ 137 138 # Test if the current data pipe exists. 139 pipes.test() 140 141 # Convert the model_from and model_to args to lists, is supplied. 142 if model_from != None: 143 model_from = [model_from] 144 if model_to != None: 145 model_to = [model_to] 146 147 # Create a list of all models. 148 models = [] 149 for model in cdp.structure.model_loop(): 150 models.append(model.num) 151 152 # Set model_from or model_to to all models if None. 153 if model_from == None: 154 model_from = models 155 if model_to == None: 156 model_to = models 157 158 # Initialise the data structure. 159 if not hasattr(cdp.structure, 'displacements'): 160 cdp.structure.displacements = Displacements() 161 162 # Loop over the starting models. 163 for i in range(len(model_from)): 164 # Assemble the atomic coordinates. 165 coord_from = [] 166 for pos in cdp.structure.atom_loop(atom_id=atom_id, model_num=model_from[i], pos_flag=True): 167 coord_from.append(pos[0]) 168 169 # Loop over the ending models. 170 for j in range(len(model_to)): 171 # Assemble the atomic coordinates. 172 coord_to = [] 173 for pos in cdp.structure.atom_loop(atom_id=atom_id, model_num=model_to[j], pos_flag=True): 174 coord_to.append(pos[0]) 175 176 # Send to the base container for the calculations. 177 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)
178 179
180 -def find_pivot(models=None, atom_id=None, init_pos=None):
181 """Superimpose a set of structural models. 182 183 @keyword models: The list of models to use. If set to None, then all models will be used. 184 @type models: list of int or None 185 @keyword atom_id: The molecule, residue, and atom identifier string. This matches the spin ID string format. 186 @type atom_id: str or None 187 @keyword init_pos: The starting pivot position for the pivot point optimisation. 188 @type init_pos: list of float or numpy rank-1, 3D array 189 """ 190 191 # Test if the current data pipe exists. 192 pipes.test() 193 194 # Initialised the starting position if needed. 195 if init_pos == None: 196 init_pos = zeros(3, float64) 197 init_pos = array(init_pos) 198 199 # Validate the models. 200 cdp.structure.validate_models() 201 202 # Create a list of all models. 203 if models == None: 204 models = [] 205 for model in cdp.structure.model_loop(): 206 models.append(model.num) 207 208 # Assemble the atomic coordinates of all models. 209 coord = [] 210 for model in models: 211 coord.append([]) 212 for pos in cdp.structure.atom_loop(atom_id=atom_id, model_num=model, pos_flag=True): 213 coord[-1].append(pos[0]) 214 coord[-1] = array(coord[-1]) 215 coord = array(coord) 216 217 # The target function. 218 finder = Pivot_finder(models, coord) 219 results = generic_minimise(func=finder.func, x0=init_pos, min_algor='simplex', print_flag=1) 220 221 # No result. 222 if results == None: 223 return 224 225 # Store the data. 226 cdp.structure.pivot = results 227 228 # Print out. 229 print("Motional pivot found at: %s" % results)
230 231
232 -def get_pos(spin_id=None, str_id=None, ave_pos=False):
233 """Load the spins from the structural object into the relax data store. 234 235 @keyword spin_id: The molecule, residue, and spin identifier string. 236 @type spin_id: str 237 @keyword str_id: The structure identifier. This can be the file name, model number, or structure number. 238 @type str_id: int or str 239 @keyword ave_pos: A flag specifying if the average atom position or the atom position from all loaded structures is loaded into the SpinContainer. 240 @type ave_pos: bool 241 """ 242 243 # Test if the current data pipe exists. 244 pipes.test() 245 246 # Test if the structure exists. 247 if not hasattr(cdp, 'structure') or not cdp.structure.num_models() or not cdp.structure.num_molecules(): 248 raise RelaxNoPdbError 249 250 # Loop over all atoms of the spin_id selection. 251 model_index = -1 252 last_model = None 253 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): 254 # Update the model info. 255 if last_model != model_num: 256 model_index = model_index + 1 257 last_model = model_num 258 259 # Remove the '+' regular expression character from the mol, res, and spin names! 260 if mol_name and search('\+', mol_name): 261 mol_name = replace(mol_name, '+', '') 262 if res_name and search('\+', res_name): 263 res_name = replace(res_name, '+', '') 264 if atom_name and search('\+', atom_name): 265 atom_name = replace(atom_name, '+', '') 266 267 # 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. 268 id = generate_spin_id(mol_name=mol_name, res_num=res_num, res_name=None, spin_name=atom_name) 269 270 # Get the spin container. 271 spin_cont = return_spin(id) 272 273 # Skip the spin if it doesn't exist. 274 if spin_cont == None: 275 continue 276 277 # Add the position vector to the spin container. 278 if ave_pos: 279 spin_cont.pos = pos 280 else: 281 if not hasattr(spin_cont, 'pos'): 282 spin_cont.pos = [] 283 spin_cont.pos.append(pos) 284 285 # Update pseudo-atoms. 286 for spin in spin_loop(): 287 if hasattr(spin, 'members'): 288 # Get the spin positions. 289 positions = [] 290 for atom in spin.members: 291 # Get the spin container. 292 subspin = return_spin(atom) 293 294 # Test that the spin exists. 295 if subspin == None: 296 raise RelaxNoSpinError(atom) 297 298 # Test the position. 299 if not hasattr(subspin, 'pos') or not subspin.pos: 300 raise RelaxError("Positional information is not available for the atom '%s'." % atom) 301 302 # Alias the position. 303 pos = subspin.pos 304 305 # Convert to a list of lists if not already. 306 multi_model = True 307 if type(pos[0]) in [float, float64]: 308 multi_model = False 309 pos = [pos] 310 311 # Store the position. 312 positions.append([]) 313 for i in range(len(pos)): 314 positions[-1].append(pos[i].tolist()) 315 316 # The averaging. 317 if spin.averaging == 'linear': 318 # Average pos. 319 ave = linear_ave(positions) 320 321 # Convert to the correct structure. 322 if multi_model: 323 spin.pos = ave 324 else: 325 spin.pos = ave[0]
326 327
328 -def load_spins(spin_id=None, str_id=None, ave_pos=False):
329 """Load the spins from the structural object into the relax data store. 330 331 @keyword spin_id: The molecule, residue, and spin identifier string. 332 @type spin_id: str 333 @keyword str_id: The structure identifier. This can be the file name, model number, or structure number. 334 @type str_id: int or str 335 @keyword ave_pos: A flag specifying if the average atom position or the atom position from all loaded structures is loaded into the SpinContainer. 336 @type ave_pos: bool 337 """ 338 339 # Test if the current data pipe exists. 340 pipes.test() 341 342 # Test if the structure exists. 343 if not hasattr(cdp, 'structure') or not cdp.structure.num_models() or not cdp.structure.num_molecules(): 344 raise RelaxNoPdbError 345 346 # Print out. 347 print("Adding the following spins to the relax data store.\n") 348 349 # Init the data for printing out. 350 mol_names = [] 351 res_nums = [] 352 res_names = [] 353 spin_nums = [] 354 spin_names = [] 355 356 # Initialise data for the atom loop. 357 model_index = -1 358 last_model = 1000000 359 360 # Loop over all atoms of the spin_id selection. 361 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): 362 # Update the model info. 363 if last_model != model_num: 364 model_index = model_index + 1 365 last_model = model_num 366 367 # Remove the '+' regular expression character from the mol, res, and spin names! 368 if mol_name and search('\+', mol_name): 369 mol_name = replace(mol_name, '+', '') 370 if res_name and search('\+', res_name): 371 res_name = replace(res_name, '+', '') 372 if atom_name and search('\+', atom_name): 373 atom_name = replace(atom_name, '+', '') 374 375 # Generate a spin ID for the current atom. 376 id = generate_spin_id(mol_name=mol_name, res_num=res_num, res_name=res_name, spin_num=atom_num, spin_name=atom_name) 377 378 # Create the spin. 379 try: 380 spin_cont = create_spin(mol_name=mol_name, res_num=res_num, res_name=res_name, spin_num=atom_num, spin_name=atom_name) 381 except RelaxError: 382 # Throw a warning if still on the first model. 383 if not model_index: 384 warn(RelaxWarning("The spin '%s' already exists." % id)) 385 continue 386 387 # Otherwise, get the spin container. 388 spin_cont = return_spin(id) 389 390 # Append all the spin ID info for the first model for printing later. 391 if model_index == 0: 392 mol_names.append(mol_name) 393 res_nums.append(res_num) 394 res_names.append(res_name) 395 spin_nums.append(atom_num) 396 spin_names.append(atom_name) 397 398 # Convert the position vector to a numpy array. 399 pos = array(pos, float64) 400 401 # Average position vector (already averaged across models in the atom_loop). 402 if ave_pos: 403 spin_cont.pos = pos 404 405 # All positions. 406 else: 407 # Initialise. 408 if not hasattr(spin_cont, 'pos'): 409 spin_cont.pos = [] 410 411 # Add the current model's position. 412 spin_cont.pos.append(pos) 413 414 # Add the element. 415 if not model_index: 416 spin_cont.element = element 417 418 # Catch no data. 419 if len(mol_names) == 0: 420 warn(RelaxWarning("No spins matching the '%s' ID string could be found." % spin_id)) 421 return 422 423 # Print out. 424 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)
425 426
427 -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):
428 """The PDB loading function. 429 430 Parsers 431 ======= 432 433 A number of parsers are available for reading PDB files. These include: 434 435 - 'scientific', the Scientific Python PDB parser. 436 - 'internal', a low quality yet fast PDB parser built into relax. 437 438 439 @keyword file: The name of the PDB file to read. 440 @type file: str 441 @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. 442 @type dir: str or None 443 @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. 444 @type read_mol: None, int, or list of int 445 @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. 446 @type set_mol_name: None, str, or list of str 447 @keyword read_model: The PDB model to extract from the file. If set to None, then all models will be loaded. 448 @type read_model: None, int, or list of int 449 @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. 450 @type set_model_num: None, int, or list of int 451 @keyword parser: The parser to be used to read the PDB file. 452 @type parser: str 453 @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. 454 @type fail: bool 455 @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. 456 @type verbosity: int 457 @raise RelaxFileError: If the fail flag is set, then a RelaxError is raised if the PDB file does not exist. 458 """ 459 460 # Test if the current data pipe exists. 461 pipes.test() 462 463 # The file path. 464 file_path = get_file_path(file, dir) 465 466 # Try adding file extensions to the end of the file path, if the file can't be found. 467 file_path_orig = file_path 468 if not access(file_path, F_OK): 469 # List of possible extensions. 470 for ext in ['.pdb', '.gz', '.pdb.gz', '.bz2', '.pdb.bz2']: 471 # Add the extension if the file can be found. 472 if access(file_path+ext, F_OK): 473 file_path = file_path + ext 474 475 # Test if the file exists. 476 if not access(file_path, F_OK): 477 if fail: 478 raise RelaxFileError('PDB', file_path_orig) 479 else: 480 warn(RelaxNoPDBFileWarning(file_path)) 481 return 482 483 # Check that the parser is the same as the currently loaded PDB files. 484 if hasattr(cdp, 'structure') and cdp.structure.id != parser: 485 raise RelaxError("The " + repr(parser) + " parser does not match the " + repr(cdp.structure.id) + " parser of the PDB loaded into the current pipe.") 486 487 # Place the parser specific structural object into the relax data store. 488 if not hasattr(cdp, 'structure'): 489 if parser == 'scientific': 490 cdp.structure = Scientific_data() 491 elif parser == 'internal': 492 cdp.structure = Internal() 493 494 # Load the structures. 495 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) 496 497 # Load into Molmol (if running). 498 molmol.molmol_obj.open_pdb()
499 500
501 -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):
502 """The XYZ loading function. 503 504 505 @keyword file: The name of the XYZ file to read. 506 @type file: str 507 @keyword dir: The directory where the XYZ file is located. If set to None, then the 508 file will be searched for in the current directory. 509 @type dir: str or None 510 @keyword read_mol: The molecule(s) to read from the file, independent of model. 511 If set to None, then all molecules will be loaded. 512 @type read_mol: None, int, or list of int 513 @keyword set_mol_name: Set the names of the molecules which are loaded. If set to None, then 514 the molecules will be automatically labelled based on the file name or 515 other information. 516 @type set_mol_name: None, str, or list of str 517 @keyword read_model: The XYZ model to extract from the file. If set to None, then all models 518 will be loaded. 519 @type read_model: None, int, or list of int 520 @keyword set_model_num: Set the model number of the loaded molecule. If set to None, then the 521 XYZ model numbers will be preserved, if they exist. 522 @type set_model_num: None, int, or list of int 523 @keyword fail: A flag which, if True, will cause a RelaxError to be raised if the XYZ 524 file does not exist. If False, then a RelaxWarning will be trown 525 instead. 526 @type fail: bool 527 @keyword verbosity: The amount of information to print to screen. Zero corresponds to 528 minimal output while higher values increase the amount of output. The 529 default value is 1. 530 @type verbosity: int 531 @raise RelaxFileError: If the fail flag is set, then a RelaxError is raised if the XYZ file 532 does not exist. 533 """ 534 535 # Test if the current data pipe exists. 536 pipes.test() 537 538 # The file path. 539 file_path = get_file_path(file, dir) 540 541 # Try adding '.xyz' to the end of the file path, if the file can't be found. 542 if not access(file_path, F_OK): 543 file_path_orig = file_path 544 file_path = file_path + '.xyz' 545 546 # Test if the file exists. 547 if not access(file_path, F_OK): 548 if fail: 549 raise RelaxFileError('XYZ', file_path_orig) 550 else: 551 warn(RelaxNoPDBFileWarning(file_path)) 552 return 553 554 # Place the structural object into the relax data store. 555 if not hasattr(cdp, 'structure'): 556 cdp.structure = Internal() 557 558 # Load the structures. 559 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)
560 561
562 -def rotate(R=None, origin=None, model=None, atom_id=None):
563 """Rotate the structural data about the origin by the specified forwards rotation. 564 565 @keyword R: The forwards rotation matrix. 566 @type R: numpy 3D, rank-2 array or a 3x3 list of floats 567 @keyword origin: The origin of the rotation. If not supplied, the origin will be set to [0, 0, 0]. 568 @type origin: numpy 3D, rank-1 array or list of len 3 or None 569 @keyword model: The model to rotate. If None, all models will be rotated. 570 @type model: int 571 @keyword atom_id: The molecule, residue, and atom identifier string. Only atoms matching this selection will be used. 572 @type atom_id: str or None 573 """ 574 575 # Test if the current data pipe exists. 576 pipes.test() 577 578 # Test if the structure exists. 579 if not hasattr(cdp, 'structure') or not cdp.structure.num_models() or not cdp.structure.num_molecules(): 580 raise RelaxNoPdbError 581 582 # Set the origin if not supplied. 583 if origin == None: 584 origin = [0, 0, 0] 585 586 # Convert the args to numpy float data structures. 587 R = array(R, float64) 588 origin = array(origin, float64) 589 590 # Call the specific code. 591 cdp.structure.rotate(R=R, origin=origin, model=model, atom_id=atom_id)
592 593
594 -def set_vector(spin=None, xh_vect=None):
595 """Place the XH unit vector into the spin container object. 596 597 @keyword spin: The spin container object. 598 @type spin: SpinContainer instance 599 @keyword xh_vect: The unit vector parallel to the XH bond. 600 @type xh_vect: array of len 3 601 """ 602 603 # Place the XH unit vector into the container. 604 spin.xh_vect = xh_vect
605 606
607 -def superimpose(models=None, method='fit to mean', atom_id=None, centroid=None):
608 """Superimpose a set of structural models. 609 610 @keyword models: The list of models to superimpose. If set to None, then all models will be used. 611 @type models: list of int or None 612 @keyword method: The superimposition method. It must be one of 'fit to mean' or 'fit to first'. 613 @type method: str 614 @keyword atom_id: The molecule, residue, and atom identifier string. This matches the spin ID string format. 615 @type atom_id: str or None 616 @keyword centroid: An alternative position of the centroid to allow for different superpositions, for example of pivot point motions. 617 @type centroid: list of float or numpy rank-1, 3D array 618 """ 619 620 # Check the method. 621 allowed = ['fit to mean', 'fit to first'] 622 if method not in allowed: 623 raise RelaxError("The superimposition method '%s' is unknown. It must be one of %s." % (method, allowed)) 624 625 # Test if the current data pipe exists. 626 pipes.test() 627 628 # Validate the models. 629 cdp.structure.validate_models() 630 631 # Create a list of all models. 632 if models == None: 633 models = [] 634 for model in cdp.structure.model_loop(): 635 models.append(model.num) 636 637 # Assemble the atomic coordinates of all models. 638 coord = [] 639 for model in models: 640 coord.append([]) 641 for pos in cdp.structure.atom_loop(atom_id=atom_id, model_num=model, pos_flag=True): 642 coord[-1].append(pos[0]) 643 coord[-1] = array(coord[-1]) 644 645 # The different algorithms. 646 if method == 'fit to mean': 647 T, R, pivot = fit_to_mean(models=models, coord=coord, centroid=centroid) 648 elif method == 'fit to first': 649 T, R, pivot = fit_to_first(models=models, coord=coord, centroid=centroid) 650 651 652 # Update to the new coordinates. 653 for i in range(len(models)): 654 # Translate the molecule first (the rotational pivot is defined in the first model). 655 translate(T=T[i], model=models[i]) 656 657 # Rotate the molecule. 658 rotate(R=R[i], origin=pivot[i], model=models[i])
659 660
661 -def translate(T=None, model=None, atom_id=None):
662 """Shift the structural data by the specified translation vector. 663 664 @keyword T: The translation vector. 665 @type T: numpy rank-1, 3D array or list of float 666 @keyword model: The model to translate. If None, all models will be rotated. 667 @type model: int or None 668 @keyword atom_id: The molecule, residue, and atom identifier string. Only atoms matching this selection will be used. 669 @type atom_id: str or None 670 """ 671 672 # Test if the current data pipe exists. 673 pipes.test() 674 675 # Test if the structure exists. 676 if not hasattr(cdp, 'structure') or not cdp.structure.num_models() or not cdp.structure.num_molecules(): 677 raise RelaxNoPdbError 678 679 # Convert the args to numpy float data structures. 680 T = array(T, float64) 681 682 # Call the specific code. 683 cdp.structure.translate(T=T, model=model, atom_id=atom_id)
684 685
686 -def vectors(attached=None, spin_id=None, model=None, verbosity=1, ave=True, unit=True):
687 """Extract the bond vectors from the loaded structures and store them in the spin container. 688 689 @keyword attached: The name of the atom attached to the spin, as given in the structural file. Regular expression can be used, for example 'H*'. This uses relax rather than Python regular expression (i.e. shell like syntax). 690 @type attached: str 691 @keyword spin_id: The spin identifier string. 692 @type spin_id: str 693 @keyword model: The model to extract the vector from. If None, all vectors will be extracted. 694 @type model: str 695 @keyword verbosity: The higher the value, the more information is printed to screen. 696 @type verbosity: int 697 @keyword ave: A flag which if True will cause the average of all vectors to be extracted. 698 @type ave: bool 699 @keyword unit: A flag which if True will cause the function to calculate the unit vectors. 700 @type unit: bool 701 """ 702 703 # Test if the current data pipe exists. 704 pipes.test() 705 706 # Test if the PDB file has been loaded. 707 if not hasattr(cdp, 'structure'): 708 raise RelaxNoPdbError 709 710 # Test if sequence data is loaded. 711 if not exists_mol_res_spin_data(): 712 raise RelaxNoSequenceError 713 714 # Print out. 715 if verbosity: 716 # Number of models. 717 num_models = cdp.structure.num_models() 718 719 # Multiple models loaded. 720 if num_models > 1: 721 if model: 722 print(("Extracting vectors for model '%s'." % model)) 723 else: 724 print(("Extracting vectors for all %s models." % num_models)) 725 if ave: 726 print("Averaging all vectors.") 727 728 # Single model loaded. 729 else: 730 print("Extracting vectors from the single model.") 731 732 # Unit vectors. 733 if unit: 734 print("Calculating the unit vectors.") 735 736 # Determine if the attached atom is a proton. 737 proton = False 738 if relax_re.search('.*H.*', attached) or relax_re.search(attached, 'H'): 739 proton = True 740 if verbosity: 741 if proton: 742 print("The attached atom is a proton.") 743 else: 744 print("The attached atom is not a proton.") 745 print('') 746 747 # Set the variable name in which the vectors will be stored. 748 if proton: 749 object_name = 'xh_vect' 750 else: 751 object_name = 'bond_vect' 752 753 # Loop over the spins. 754 no_vectors = True 755 for spin, mol_name, res_num, res_name in spin_loop(selection=spin_id, full_info=True): 756 # Skip deselected spins. 757 if not spin.select: 758 continue 759 760 # 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. 761 id = generate_spin_id(res_num=res_num, res_name=None, spin_name=spin.name, spin_num=spin.num) 762 763 # Test that the spin number or name are set (one or both are essential for the identification of the atom). 764 if spin.num == None and spin.name == None: 765 warn(RelaxWarning("Either the spin number or name must be set for the spin " + repr(id) + " to identify the corresponding atom in the molecule.")) 766 continue 767 768 # The bond vector already exists. 769 if hasattr(spin, object_name): 770 obj = getattr(spin, object_name) 771 if obj != None: 772 warn(RelaxWarning("The bond vector for the spin " + repr(id) + " already exists.")) 773 continue 774 775 # Get the bond info. 776 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) 777 778 # No attached atom. 779 if not bond_vectors: 780 # Warning messages. 781 if warnings: 782 warn(RelaxWarning(warnings + " (atom ID " + repr(id) + ").")) 783 784 # Skip the spin. 785 continue 786 787 # Set the attached atom name. 788 if not hasattr(spin, 'attached_atom'): 789 spin.attached_atom = attached_name 790 elif spin.attached_atom != attached_name: 791 raise RelaxError("The " + repr(spin.attached_atom) + " atom already attached to the spin does not match the attached atom " + repr(attached_name) + ".") 792 793 # Initialise the average vector. 794 if ave: 795 ave_vector = zeros(3, float64) 796 797 # Loop over the individual vectors. 798 for i in xrange(len(bond_vectors)): 799 # Unit vector. 800 if unit: 801 # Normalisation factor. 802 norm_factor = norm(bond_vectors[i]) 803 804 # Test for zero length. 805 if norm_factor == 0.0: 806 warn(RelaxZeroVectorWarning(id)) 807 808 # Calculate the normalised vector. 809 else: 810 bond_vectors[i] = bond_vectors[i] / norm_factor 811 812 # Sum the vectors. 813 if ave: 814 ave_vector = ave_vector + bond_vectors[i] 815 816 # Average. 817 if ave: 818 vector = ave_vector / float(len(bond_vectors)) 819 else: 820 vector = bond_vectors 821 822 # Convert to a single vector if needed. 823 if len(vector) == 1: 824 vector = vector[0] 825 826 # Set the vector. 827 setattr(spin, object_name, vector) 828 829 # We have a vector! 830 no_vectors = False 831 832 # Print out of modified spins. 833 if verbosity: 834 # The number of vectors. 835 num = len(bond_vectors) 836 plural = 's' 837 if num == 1: 838 plural = '' 839 840 if spin.name: 841 print("Extracted %s %s-%s vector%s for the spin '%s'." % (num, spin.name, attached_name, plural, id)) 842 else: 843 print("Extracted %s %s-%s vector%s for the spin '%s'." % (num, spin.num, attached_name, plural, id)) 844 845 # Right, catch the problem of missing vectors to prevent massive user confusion! 846 if no_vectors: 847 raise RelaxError("No vectors could be extracted.")
848 849
850 -def write_pdb(file=None, dir=None, model_num=None, compress_type=0, force=False):
851 """The PDB writing function. 852 853 @keyword file: The name of the PDB file to write. 854 @type file: str 855 @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. 856 @type dir: str or None 857 @keyword model_num: The model to place into the PDB file. If not supplied, then all models will be placed into the file. 858 @type model_num: None or int 859 @keyword compress_type: The compression type. The integer values correspond to the compression type: 0, no compression; 1, Bzip2 compression; 2, Gzip compression. 860 @type compress_type: int 861 @keyword force: The force flag which if True will cause the file to be overwritten. 862 @type force: bool 863 """ 864 865 # Test if the current data pipe exists. 866 pipes.test() 867 868 # Check if the structural object exists. 869 if not hasattr(cdp, 'structure'): 870 raise RelaxError("No structural data is present in the current data pipe.") 871 872 # Check if the structural object is writable. 873 if cdp.structure.id in ['scientific']: 874 raise RelaxError("The structures from the " + cdp.structure.id + " parser are not writable.") 875 876 # The file path. 877 file_path = get_file_path(file, dir) 878 879 # Add '.pdb' to the end of the file path if it isn't there yet. 880 if not search(".pdb$", file_path): 881 file_path = file_path + '.pdb' 882 883 # Open the file for writing. 884 file = open_write_file(file_path, compress_type=compress_type, force=force) 885 886 # Write the structures. 887 cdp.structure.write_pdb(file, model_num=model_num)
888