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

Source Code for Module generic_fns.structure.internal

   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  # Module docstring. 
  23  """Module containing the internal relax structural object.""" 
  24   
  25  # Python module imports. 
  26  from numpy import array, dot, float64, linalg, zeros 
  27  import os 
  28  from os import F_OK, access 
  29  from re import search 
  30  from string import digits 
  31  from warnings import warn 
  32   
  33  # relax module imports. 
  34  from data.relax_xml import fill_object_contents, xml_to_object 
  35  from generic_fns import pipes, relax_re 
  36  from generic_fns.mol_res_spin import spin_loop 
  37  from generic_fns.mol_res_spin import Selection 
  38  from generic_fns.structure.api_base import Base_struct_API, ModelList, Displacements 
  39  from relax_errors import RelaxError, RelaxNoneIntError, RelaxNoPdbError 
  40  from relax_io import file_root, open_read_file 
  41  from relax_warnings import RelaxWarning 
  42   
  43   
  44   
45 -class Internal(Base_struct_API):
46 """The internal relax structural data object. 47 48 The structural data object for this class is a container possessing a number of different arrays 49 corresponding to different structural information. These objects are described in the 50 structural container docstring. 51 """ 52 53 # Identification string. 54 id = 'internal' 55 56
57 - def _bonded_atom(self, attached_atom, index, mol):
58 """Find the atom named attached_atom directly bonded to the atom located at the index. 59 60 @param attached_atom: The name of the attached atom to return. 61 @type attached_atom: str 62 @param index: The index of the atom which the attached atom is attached to. 63 @type index: int 64 @param mol: The molecule container. 65 @type mol: MolContainer instance 66 @return: A tuple of information about the bonded atom. 67 @rtype: tuple consisting of the atom number (int), atom name (str), element 68 name (str), and atomic position (Numeric array of len 3) 69 """ 70 71 # Init. 72 bonded_found = False 73 74 # No bonded atoms, so determine the connectivities. 75 if not mol.bonded[index]: 76 # Determine the molecule type if needed. 77 if not hasattr(mol, 'type'): 78 self._mol_type(mol) 79 80 # Protein. 81 if mol.type == 'protein': 82 self._protein_connect(mol) 83 84 # Find everything within 2 Angstroms and say they are bonded. 85 else: 86 self._find_bonded_atoms(index, mol, radius=2) 87 88 # Loop over the bonded atoms. 89 matching_list = [] 90 for bonded_index in mol.bonded[index]: 91 if relax_re.search(mol.atom_name[bonded_index], attached_atom): 92 matching_list.append(bonded_index) 93 num_attached = len(matching_list) 94 95 # Problem. 96 if num_attached > 1: 97 # Get the atom names. 98 matching_names = [] 99 for i in matching_list: 100 matching_names.append(mol.atom_name[i]) 101 102 # Return nothing but a warning. 103 return None, None, None, None, None, 'More than one attached atom found: ' + repr(matching_names) 104 105 # No attached atoms. 106 if num_attached == 0: 107 if relax_re.search('@*', attached_atom): 108 matching_list = [] 109 bonded_num=[] 110 bonded_name=[] 111 element=[] 112 pos=[] 113 for spin, mol_name, res_num, res_name in spin_loop(selection=attached_atom, full_info=True): 114 bonded_num.append(spin.num) 115 bonded_name.append(spin.name) 116 element.append(spin.element) 117 pos.append(spin.pos) 118 if len(bonded_num) == 1: 119 return bonded_num[0], bonded_name[0], element[0], pos[0], attached_atom, None 120 elif len(bonded_num) > 1: 121 # Return nothing but a warning. 122 return None, None, None, None, None, 'More than one attached atom found: ' + repr(matching_names) 123 elif len(bonded_num) > 1: 124 # Return nothing but a warning. 125 return None, None, None, None, None, "No attached atom could be found" 126 else: 127 return None, None, None, None, None, "No attached atom could be found" 128 129 # The bonded atom info. 130 index = matching_list[0] 131 bonded_num = mol.atom_num[index] 132 bonded_name = mol.atom_name[index] 133 element = mol.element[index] 134 pos = [mol.x[index], mol.y[index], mol.z[index]] 135 attached_name = mol.atom_name[index] 136 137 # Return the information. 138 return bonded_num, bonded_name, element, pos, attached_name, None
139 140
141 - def _find_bonded_atoms(self, index, mol, radius=1.2):
142 """Find all atoms within a sphere and say that they are attached to the central atom. 143 144 The found atoms will be added to the 'bonded' data structure. 145 146 147 @param index: The index of the central atom. 148 @type index: int 149 @param mol: The molecule container. 150 @type mol: MolContainer instance 151 """ 152 153 # Central atom info. 154 centre = array([mol.x[index], mol.y[index], mol.z[index]], float64) 155 156 # Atom loop. 157 dist_list = [] 158 connect_list = {} 159 element_list = {} 160 for i in range(len(mol.atom_num)): 161 # Skip proton to proton bonds! 162 if mol.element[index] == 'H' and mol.element[i] == 'H': 163 continue 164 165 # The atom's position. 166 pos = array([mol.x[i], mol.y[i], mol.z[i]], float64) 167 168 # The distance from the centre. 169 dist = linalg.norm(centre-pos) 170 171 # The atom is within the radius. 172 if dist < radius: 173 # Store the distance. 174 dist_list.append(dist) 175 176 # Store the atom index. 177 connect_list[dist] = i 178 179 # Store the element type. 180 element_list[dist] = mol.element[i] 181 182 # The maximum number of allowed covalent bonds. 183 max_conn = 1000 # Ridiculous default! 184 if mol.element[index] == 'H': 185 max_conn = 1 186 elif mol.element[index] == 'O': 187 max_conn = 2 188 elif mol.element[index] == 'N': 189 max_conn = 3 190 elif mol.element[index] == 'C': 191 max_conn = 4 192 193 # Sort. 194 dist_list.sort() 195 196 # Loop over the max number of connections (or the number of connected atoms, if less). 197 for i in range(min(max_conn, len(dist_list))): 198 mol.atom_connect(index, connect_list[dist_list[i]])
199 200
201 - def _get_chemical_name(self, hetID):
202 """Return the chemical name corresponding to the given residue ID. 203 204 The following names are currently returned:: 205 ________________________________________________ 206 | | | 207 | hetID | Chemical name | 208 |________|_____________________________________| 209 | | | 210 | TNS | Tensor | 211 | COM | Centre of mass | 212 | AXS | Tensor axes | 213 | SIM | Monte Carlo simulation tensor axes | 214 | PIV | Pivot point | 215 | CON | Cone object | 216 | AVE | Average vector | 217 |________|_____________________________________| 218 219 For any other residues, no description is returned. 220 221 @param hetID: The residue ID. 222 @type hetID: str 223 @return: The chemical name. 224 @rtype: str or None 225 """ 226 227 # Tensor. 228 if hetID == 'TNS': 229 return 'Tensor' 230 231 # Centre of mass. 232 if hetID == 'COM': 233 return 'Centre of mass' 234 235 # Tensor axes. 236 if hetID == 'AXS': 237 return 'Tensor axes' 238 239 # Monte Carlo simulation tensor axes. 240 if hetID == 'SIM': 241 return 'Monte Carlo simulation tensor axes' 242 243 # Pivot point. 244 if hetID == 'PIV': 245 return 'Pivot point' 246 247 # Cone object. 248 if hetID == 'CON': 249 return 'Cone' 250 251 # Average vector. 252 if hetID == 'AVE': 253 return 'Average vector'
254 255
256 - def _parse_models_pdb(self, file_path):
257 """Generator function for looping over the models in the PDB file. 258 259 @param file_path: The full path of the PDB file. 260 @type file_path: str 261 @return: The model number and all the records for that model. 262 @rtype: tuple of int and array of str 263 """ 264 265 # Open the file. 266 file = open_read_file(file_path) 267 lines = file.readlines() 268 file.close() 269 270 # Check for empty files. 271 if lines == []: 272 raise RelaxError("The PDB file is empty.") 273 274 # Init. 275 model = None 276 records = [] 277 278 # Loop over the data. 279 for i in range(len(lines)): 280 # A new model record. 281 if search('^MODEL', lines[i]): 282 try: 283 model = int(lines[i].split()[1]) 284 except: 285 raise RelaxError("The MODEL record " + repr(lines[i]) + " is corrupt, cannot read the PDB file.") 286 287 # Skip all records prior to the first ATOM or HETATM record. 288 if not (search('^ATOM', lines[i]) or search('^HETATM', lines[i])) and not len(records): 289 continue 290 291 # End of the model. 292 if search('^ENDMDL', lines[i]): 293 # Yield the info. 294 yield model, records 295 296 # Reset the records. 297 records = [] 298 299 # Skip the rest of this loop. 300 continue 301 302 # Append the line as a record of the model. 303 records.append(lines[i]) 304 305 # If records is not empty then there are no models, so yield the lot. 306 if len(records): 307 yield model, records
308 309
310 - def _parse_models_xyz(self, file_path):
311 """Generator function for looping over the models in the XYZ file. 312 313 @param file_path: The full path of the XYZ file. 314 @type file_path: str 315 @return: The model number and all the records for that model. 316 @rtype: tuple of int and array of str 317 """ 318 319 # Open the file. 320 file = open_read_file(file_path) 321 lines = file.readlines() 322 file.close() 323 324 # Check for empty files. 325 if lines == []: 326 raise RelaxError("The XYZ file is empty.") 327 328 # Init. 329 total_atom = 0 330 model = 0 331 records = [] 332 333 # Loop over the data. 334 for i in range(len(lines)): 335 num=0 336 word = lines[i].split() 337 # Find the total atom number and the first model. 338 if (i==0) and (len(word)==1): 339 try: 340 total_atom = int(word[0]) 341 num = 1 342 except: 343 raise RelaxError("The MODEL record " + repr(lines[i]) + " is corrupt, cannot read the XYZ file.") 344 345 # End of the model. 346 if (len(records) == total_atom): 347 # Yield the info 348 yield records 349 350 # Reset the records. 351 records = [] 352 353 # Skip all records prior to atom coordinates record. 354 if (len(word) != 4): 355 continue 356 357 # Append the line as a record of the model. 358 records.append(lines[i]) 359 360 # If records is not empty then there are no models, so yield the lot. 361 if len(records): 362 yield records
363 364
365 - def _parse_mols(self, records):
366 """Generator function for looping over the molecules in the PDB records of a model. 367 368 @param records: The list of PDB records for the model, or if no models exist the entire 369 PDB file. 370 @type records: list of str 371 @return: The molecule number and all the records for that molecule. 372 @rtype: tuple of int and list of str 373 """ 374 375 # Check for empty records. 376 if records == []: 377 raise RelaxError("There are no PDB records for this model.") 378 379 # Init. 380 mol_num = 1 381 mol_records = [] 382 end = False 383 384 # Loop over the data. 385 for i in range(len(records)): 386 # A PDB termination record. 387 if search('^END', records[i]): 388 break 389 390 # A model termination record. 391 if search('^ENDMDL', records[i]): 392 end = True 393 394 # A molecule termination record with no trailing HETATM. 395 elif i < len(records)-1 and search('^TER', records[i]) and not search('^HETATM', records[i+1]): 396 end = True 397 398 # A HETATM followed by an ATOM record. 399 elif i < len(records)-1 and search('^HETATM', records[i]) and search('^ATOM', records[i+1]): 400 end = True 401 402 # End. 403 if end: 404 # Yield the info. 405 yield mol_num, mol_records 406 407 # Reset the records. 408 mol_records = [] 409 410 # Increment the molecule number. 411 mol_num = mol_num + 1 412 413 # Reset the flag. 414 end = False 415 416 # Skip the rest of this loop. 417 continue 418 419 # Append the line as a record of the molecule. 420 mol_records.append(records[i]) 421 422 # If records is not empty then there is only a single molecule, so yield the lot. 423 if len(mol_records): 424 yield mol_num, mol_records
425 426
427 - def _validate_data_arrays(self, struct):
428 """Check the validity of the data arrays in the given structure object. 429 430 @param struct: The structural object. 431 @type struct: Structure_container instance 432 """ 433 434 # The number of atoms. 435 num = len(struct.atom_name) 436 437 # Check the other lengths. 438 if len(struct.bonded) != num and len(struct.chain_id) != num and len(struct.element) != num and len(struct.pdb_record) != num and len(struct.res_name) != num and len(struct.res_num) != num and len(struct.seg_id) != num and len(struct.x) != num and len(struct.y) != num and len(struct.z) != num: 439 raise RelaxError("The structural data is invalid.")
440 441
442 - def _mol_type(self, mol):
443 """Determine the type of molecule. 444 445 @param mol: The molecule data container. 446 @type mol: MolContainer instance 447 """ 448 449 # Amino acids. 450 aa = ['ALA', 'ARG', 'ASN', 'ASP', 'CYS', 'GLU', 'GLN', 'GLY', 'HIS', 'ILE', 'LEU', 'LYS', 'MET', 'PHE', 'PRO', 'SER', 'THR', 'TRP', 'TYR', 'VAL'] 451 452 # Set the molecule type to default to 'other'. 453 mol.type = 'other' 454 455 # Loop over the residues. 456 for res in mol.res_name: 457 # Protein. 458 if res in aa: 459 # Set the molecule type and return. 460 mol.type = 'protein' 461 return
462 463
464 - def _protein_connect(self, mol):
465 """Set up the connectivities for the protein. 466 467 @param mol: The molecule data container. 468 @type mol: MolContainer instance 469 """ 470 471 # Initialise some residue data. 472 curr_res_num = None 473 res_atoms = [] 474 475 # Loop over all atoms. 476 for i in range(len(mol.atom_num)): 477 # New residue. 478 if mol.res_num[i] != curr_res_num: 479 # Intra-residue connectivites. 480 if len(res_atoms): 481 self._protein_intra_connect(mol, res_atoms) 482 483 # Update the residue number. 484 curr_res_num = mol.res_num[i] 485 486 # Reset the residue atom index list. 487 res_atoms = [] 488 489 # Add the atom index to the list. 490 res_atoms.append(i) 491 492 # Last atom. 493 if i == len(mol.atom_num) - 1 and len(res_atoms): 494 self._protein_intra_connect(mol, res_atoms)
495 496
497 - def _protein_intra_connect(self, mol, res_atoms):
498 """Set up the connectivities for the protein. 499 500 @param mol: The molecule data container. 501 @type mol: MolContainer instance 502 @param res_atoms: The list of atom indices corresponding to the residue. 503 @type res_atoms: list of int 504 """ 505 506 # Back bond connectivity. 507 indices = { 508 'N': None, 509 'C': None, 510 'O': None, 511 'CA': None, 512 'HN': None, 513 'H': None, # Same as HN. 514 'HA': None 515 } 516 517 # Loop over all atoms to find the indices. 518 for index in res_atoms: 519 if mol.atom_name[index] in indices: 520 indices[mol.atom_name[index]] = index 521 522 # Connect the atom pairs. 523 pairs = [ 524 ['N', 'HN'], 525 ['N', 'H'], 526 ['N', 'CA'], 527 ['CA', 'HA'], 528 ['CA', 'C'], 529 ['C', 'O'] 530 ] 531 532 # Loop over the atoms pairs and connect them. 533 for pair in pairs: 534 if indices[pair[0]] != None and indices[pair[1]] != None: 535 mol.atom_connect(indices[pair[0]], indices[pair[1]])
536 537
538 - def _translate(self, data, format='str'):
539 """Convert the data into a format for writing to file. 540 541 @param data: The data to convert to the required format. 542 @type data: anything 543 @keyword format: The format to convert to. This can be 'str', 'float', or 'int'. 544 @type format: str 545 @return: The converted version of the data. 546 @rtype: str 547 """ 548 549 # Conversion to string. 550 if format == 'str': 551 # None values. 552 if data == None: 553 data = '' 554 555 # Force convert to string. 556 if not isinstance(data, str): 557 data = repr(data) 558 559 # Conversion to float. 560 if format == 'float': 561 # None values. 562 if data == None: 563 data = 0.0 564 565 # Force convert to float. 566 if not isinstance(data, float): 567 data = float(data) 568 569 # Return the converted data. 570 return data
571 572
573 - def add_atom(self, 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):
574 """Add a new atom to the structural data object. 575 576 @keyword mol_name: The name of the molecule. 577 @type mol_name: str 578 @keyword atom_name: The atom name, e.g. 'H1'. 579 @type atom_name: str or None 580 @keyword res_name: The residue name. 581 @type res_name: str or None 582 @keyword res_num: The residue number. 583 @type res_num: int or None 584 @keyword pos: The position vector of coordinates. 585 @type pos: list (length = 3) 586 @keyword element: The element symbol. 587 @type element: str or None 588 @keyword atom_num: The atom number. 589 @type atom_num: int or None 590 @keyword chain_id: The chain identifier. 591 @type chain_id: str or None 592 @keyword segment_id: The segment identifier. 593 @type segment_id: str or None 594 @keyword pdb_record: The optional PDB record name, e.g. 'ATOM' or 'HETATM'. 595 @type pdb_record: str or None 596 """ 597 598 # Test if the current data pipe exists. 599 pipes.test() 600 601 # Add a model if not present. 602 if len(self.structural_data) == 0: 603 self.add_model() 604 605 # Loop over each model. 606 for model in self.structural_data: 607 # Specific molecule. 608 mol = self.get_molecule(mol_name, model=model.num) 609 610 # Add the molecule, if it does not exist. 611 if mol == None: 612 self.add_molecule(name=mol_name) 613 mol = self.get_molecule(mol_name, model=model.num) 614 615 # Add the atom. 616 mol.atom_add(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)
617 618
619 - def add_model(self, model=None, coords_from=None):
620 """Add a new model to the store. 621 622 The new model will be constructured with the structural information from the other models currently present. The coords_from argument allows the atomic positions to be taken from a certain model. If this argument is not set, then the atomic positions from the first model will be used. 623 624 @keyword model: The number of the model to create. 625 @type model: int or None 626 @keyword coords_from: The model number to take the coordinates from. 627 @type coords_from: int or None 628 @return: The model container. 629 @rtype: ModelContainer instance 630 """ 631 632 # Check if the model currently exists. 633 if model != None: 634 for i in range(len(self.structural_data)): 635 if model == self.structural_data[i].num: 636 raise RelaxError("The model '%s' already exists." % model) 637 638 # Add a new model. 639 self.structural_data.add_item(model_num=model) 640 641 # The model to duplicate. 642 if coords_from == None: 643 coords_from = self.structural_data[0].num 644 645 # Construct the structural data for the model from the other models. 646 for mol_name, res_num, res_name, atom_num, atom_name, element, pos in self.atom_loop(model_num=coords_from, 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): 647 # Add the atom. 648 self.add_atom(self, mol_name=mol_name, atom_name=atom_name, res_name=res_name, res_num=res_num, pos=pos, element=element, atom_num=atom_num) 649 650 # Return the model. 651 return self.structural_data[-1]
652 653
654 - def add_molecule(self, name=None):
655 """Add a new molecule to the store. 656 657 @keyword name: The molecule identifier string. 658 @type name: str 659 """ 660 661 # Add a model if necessary. 662 if len(self.structural_data) == 0: 663 self.add_model() 664 665 # Loop over the models. 666 for i in range(len(self.structural_data)): 667 # Add the molecule. 668 self.structural_data[i].mol.add_item(mol_name=name, mol_cont=MolContainer())
669 670
671 - def are_bonded(self, atom_id1=None, atom_id2=None):
672 """Determine if two atoms are directly bonded to each other. 673 674 @keyword atom_id1: The molecule, residue, and atom identifier string of the first atom. 675 @type atom_id1: str 676 @keyword atom_id2: The molecule, residue, and atom identifier string of the second atom. 677 @type atom_id2: str 678 @return: True if the atoms are directly bonded. 679 @rtype: bool 680 """ 681 682 # Generate the selection objects. 683 sel_obj1 = Selection(atom_id1) 684 sel_obj2 = Selection(atom_id2) 685 686 # Build the connectivities if needed. 687 for mol in self.structural_data[0].mol: 688 for i in range(len(mol.atom_num)): 689 if not len(mol.bonded[i]): 690 self._find_bonded_atoms(i, mol, radius=2) 691 692 # Loop over the molecules. 693 for mol in self.structural_data[0].mol: 694 # Skip non-matching molecules. 695 if not sel_obj1.contains_mol(mol.mol_name): 696 continue 697 if not sel_obj2.contains_mol(mol.mol_name): 698 continue 699 700 # Find the first atom. 701 index1 = None 702 for i in range(len(mol.atom_num)): 703 # Skip a non-matching first atom. 704 if sel_obj1.contains_spin(mol.atom_num[i], mol.atom_name[i], mol.res_num[i], mol.res_name[i], mol.mol_name): 705 index1 = i 706 break 707 708 # Find the second atom. 709 index2 = None 710 for i in range(len(mol.atom_num)): 711 # Skip a non-matching first atom. 712 if sel_obj2.contains_spin(mol.atom_num[i], mol.atom_name[i], mol.res_num[i], mol.res_name[i], mol.mol_name): 713 index2 = i 714 break 715 716 # Connectivities exist. 717 if index1 < len(mol.bonded): 718 if index2 in mol.bonded[index1]: 719 return True 720 else: 721 return False
722 723
724 - def atom_loop(self, atom_id=None, str_id=None, model_num=None, model_num_flag=False, mol_name_flag=False, res_num_flag=False, res_name_flag=False, atom_num_flag=False, atom_name_flag=False, element_flag=False, pos_flag=False, ave=False):
725 """Generator function for looping over all atoms in the internal relax structural object. 726 727 @keyword atom_id: The molecule, residue, and atom identifier string. Only atoms matching this selection will be yielded. 728 @type atom_id: str 729 @keyword str_id: The structure identifier. This can be the file name, model number, or structure number. If None, then all structures will be looped over. 730 @type str_id: str, int, or None 731 @keyword model_num: Only loop over a specific model. 732 @type model_num: int or None 733 @keyword model_num_flag: A flag which if True will cause the model number to be yielded. 734 @type model_num_flag: bool 735 @keyword mol_name_flag: A flag which if True will cause the molecule name to be yielded. 736 @type mol_name_flag: bool 737 @keyword res_num_flag: A flag which if True will cause the residue number to be yielded. 738 @type res_num_flag: bool 739 @keyword res_name_flag: A flag which if True will cause the residue name to be yielded. 740 @type res_name_flag: bool 741 @keyword atom_num_flag: A flag which if True will cause the atom number to be yielded. 742 @type atom_num_flag: bool 743 @keyword atom_name_flag: A flag which if True will cause the atom name to be yielded. 744 @type atom_name_flag: bool 745 @keyword element_flag: A flag which if True will cause the element name to be yielded. 746 @type element_flag: bool 747 @keyword pos_flag: A flag which if True will cause the atomic position to be yielded. 748 @type pos_flag: bool 749 @keyword ave: A flag which if True will result in this method returning the average atom properties across all loaded structures. 750 @type ave: bool 751 @return: A tuple of atomic information, as described in the docstring. 752 @rtype: tuple consisting of optional molecule name (str), residue number (int), residue name (str), atom number (int), atom name(str), element name (str), and atomic position (array of len 3). 753 """ 754 755 # Check that the structure is loaded. 756 if not len(self.structural_data): 757 raise RelaxNoPdbError 758 759 # Generate the selection object. 760 sel_obj = None 761 if atom_id: 762 sel_obj = Selection(atom_id) 763 764 # Model loop. 765 for model in self.model_loop(model_num): 766 # Loop over the molecules. 767 for mol_index in range(len(model.mol)): 768 mol = model.mol[mol_index] 769 770 # Skip non-matching molecules. 771 if sel_obj and not sel_obj.contains_mol(mol.mol_name): 772 continue 773 774 # Loop over all atoms. 775 for i in range(len(mol.atom_name)): 776 # Skip non-matching atoms. 777 if sel_obj and not sel_obj.contains_spin(mol.atom_num[i], mol.atom_name[i], mol.res_num[i], mol.res_name[i], mol.mol_name): 778 continue 779 780 # Initialise. 781 res_num = mol.res_num[i] 782 res_name = mol.res_name[i] 783 atom_num = mol.atom_num[i] 784 atom_name = mol.atom_name[i] 785 element = mol.element[i] 786 pos = zeros(3, float64) 787 788 # The atom position. 789 if ave: 790 # Loop over the models. 791 for model in self.model_loop(): 792 # Alias. 793 mol = model.mol[mol_index] 794 795 # Some sanity checks. 796 if mol.atom_num[i] != atom_num: 797 raise RelaxError("The loaded structures do not contain the same atoms. The average structural properties can not be calculated.") 798 799 # Sum the atom positions. 800 pos = pos + array([mol.x[i], mol.y[i], mol.z[i]], float64) 801 802 # Average the position array (divide by the number of models). 803 pos = pos / len(self.structural_data) 804 else: 805 pos = array([mol.x[i], mol.y[i], mol.z[i]], float64) 806 807 # The molecule name. 808 mol_name = mol.mol_name 809 810 # Build the tuple to be yielded. 811 atomic_tuple = () 812 if model_num_flag: 813 if ave: 814 atomic_tuple = atomic_tuple + (None,) 815 else: 816 atomic_tuple = atomic_tuple + (model.num,) 817 if mol_name_flag: 818 atomic_tuple = atomic_tuple + (mol_name,) 819 if res_num_flag: 820 atomic_tuple = atomic_tuple + (res_num,) 821 if res_name_flag: 822 atomic_tuple = atomic_tuple + (res_name,) 823 if atom_num_flag: 824 atomic_tuple = atomic_tuple + (atom_num,) 825 if atom_name_flag: 826 atomic_tuple = atomic_tuple + (atom_name,) 827 if element_flag: 828 atomic_tuple = atomic_tuple + (element,) 829 if pos_flag: 830 atomic_tuple = atomic_tuple + (pos,) 831 832 # Yield the information. 833 yield atomic_tuple 834 835 # Break out of the loop if the ave flag is set, as data from only one model is used. 836 if ave: 837 break
838 839
840 - def bond_vectors(self, attached_atom=None, model_num=None, mol_name=None, res_num=None, res_name=None, spin_num=None, spin_name=None, return_name=False, return_warnings=False):
841 """Find the bond vector between the atoms of 'attached_atom' and 'atom_id'. 842 843 @keyword attached_atom: The name of the bonded atom. 844 @type attached_atom: str 845 @keyword model_num: The model of which to return the vectors from. If not supplied and multiple models exist, then vectors from all models will be returned. 846 @type model_num: None or int 847 @keyword mol_name: The name of the molecule that attached_atom belongs to. 848 @type mol_name: str 849 @keyword res_num: The number of the residue that attached_atom belongs to. 850 @type res_num: str 851 @keyword res_name: The name of the residue that attached_atom belongs to. 852 @type res_name: str 853 @keyword spin_num: The number of the spin that attached_atom is attached to. 854 @type spin_num: str 855 @keyword spin_name: The name of the spin that attached_atom is attached to. 856 @type spin_name: str 857 @keyword return_name: A flag which if True will cause the name of the attached atom to be returned together with the bond vectors. 858 @type return_name: bool 859 @keyword return_warnings: A flag which if True will cause warning messages to be returned. 860 @type return_warnings: bool 861 @return: The list of bond vectors for each model. 862 @rtype: list of numpy arrays (or a tuple if return_name or return_warnings are set) 863 """ 864 865 # Initialise some objects. 866 vectors = [] 867 attached_name = None 868 warnings = None 869 870 # Loop over the models. 871 for model in self.model_loop(model_num): 872 # Loop over the molecules. 873 for mol in model.mol: 874 # Skip non-matching molecules. 875 if mol_name and mol_name != mol.mol_name: 876 continue 877 878 # Find the atomic index of the base atom. 879 index = None 880 for i in range(len(mol.atom_name)): 881 # Residues don't match. 882 if (res_num != None and mol.res_num[i] != res_num) or (res_name != None and mol.res_name[i] != res_name): 883 continue 884 885 # Atoms don't match. 886 if (spin_num != None and mol.atom_num[i] != spin_num) or (spin_name != None and mol.atom_name[i] != spin_name): 887 continue 888 889 # Update the index and stop searching. 890 index = i 891 break 892 893 # Found the atom. 894 if index != None: 895 # Get the atom bonded to this model/molecule/residue/atom. 896 bonded_num, bonded_name, element, pos, attached_name, warnings = self._bonded_atom(attached_atom, index, mol) 897 898 # No bonded atom. 899 if (bonded_num, bonded_name, element) == (None, None, None): 900 continue 901 902 # The bond vector. 903 vector = array(pos, float64) - array([mol.x[index], mol.y[index], mol.z[index]], float64) 904 905 # Append the vector to the vectors array. 906 vectors.append(vector) 907 908 # Not found. 909 else: 910 warnings = "Cannot find the atom in the structure" 911 912 # Build the tuple to be yielded. 913 data = (vectors,) 914 if return_name: 915 data = data + (attached_name,) 916 if return_warnings: 917 data = data + (warnings,) 918 919 # Return the data. 920 return data
921 922
923 - def connect_atom(self, mol_name=None, index1=None, index2=None):
924 """Connect two atoms in the structural data object. 925 926 @keyword mol_name: The name of the molecule. 927 @type mol_name: str 928 @keyword index1: The global index of the first atom. 929 @type index1: str 930 @keyword index2: The global index of the first atom. 931 @type index2: str 932 """ 933 934 # Test if the current data pipe exists. 935 pipes.test() 936 937 # Add the molecule, if it does not exist. 938 if self.get_molecule(mol_name) == None: 939 self.add_molecule(name=mol_name) 940 941 # Loop over each model. 942 for model in self.structural_data: 943 # Specific molecule. 944 mol = self.get_molecule(mol_name) 945 946 # Add the atom. 947 mol.atom_connect(index1=index1, index2=index2)
948 949
950 - def delete(self):
951 """Delete all the structural information.""" 952 953 # Print out. 954 print("Deleting the following structural data:\n") 955 print(self.structural_data) 956 957 # Delete the structural data. 958 del self.structural_data 959 960 # Initialise the empty model list. 961 self.structural_data = ModelList()
962 963
964 - def get_molecule(self, molecule, model=None):
965 """Return the molecule. 966 967 Only one model can be specified. 968 969 970 @param molecule: The molecule name. 971 @type molecule: int or None 972 @keyword model: The model number. 973 @type model: int or None 974 @raises RelaxError: If the model is not specified and there is more than one model loaded. 975 @return: The MolContainer corresponding to the molecule name and model number. 976 @rtype: MolContainer instance or None 977 """ 978 979 # Check if the target is a single molecule. 980 if model == None and self.num_models() > 1: 981 raise RelaxError("The target molecule cannot be determined as there are %s models already present." % self.num_models()) 982 983 # Check the model argument. 984 if not isinstance(model, int) and not model == None: 985 raise RelaxNoneIntError 986 987 # No models. 988 if not len(self.structural_data): 989 return 990 991 # Loop over the models. 992 for model_cont in self.model_loop(model): 993 # Loop over the molecules. 994 for mol in model_cont.mol: 995 # Return the matching molecule. 996 if mol.mol_name == molecule: 997 return mol
998 999
1000 - def load_pdb(self, file_path, read_mol=None, set_mol_name=None, read_model=None, set_model_num=None, verbosity=False):
1001 """Method for loading structures from a PDB file. 1002 1003 @param file_path: The full path of the PDB file. 1004 @type file_path: str 1005 @keyword read_mol: The molecule(s) to read from the file, independent of model. The 1006 molecules are determined differently by the different parsers, but 1007 are numbered consecutively from 1. If set to None, then all 1008 molecules will be loaded. 1009 @type read_mol: None, int, or list of int 1010 @keyword set_mol_name: Set the names of the molecules which are loaded. If set to None, 1011 then the molecules will be automatically labelled based on the file 1012 name or other information. 1013 @type set_mol_name: None, str, or list of str 1014 @keyword read_model: The PDB model to extract from the file. If set to None, then all 1015 models will be loaded. 1016 @type read_model: None, int, or list of int 1017 @keyword set_model_num: Set the model number of the loaded molecule. If set to None, then 1018 the PDB model numbers will be preserved, if they exist. 1019 @type set_model_num: None, int, or list of int 1020 @keyword verbosity: A flag which if True will cause messages to be printed. 1021 @type verbosity: bool 1022 @return: The status of the loading of the PDB file. 1023 @rtype: bool 1024 """ 1025 1026 # Initial printout. 1027 if verbosity: 1028 print("\nInternal relax PDB parser.") 1029 1030 # Test if the file exists. 1031 if not access(file_path, F_OK): 1032 # Exit indicating failure. 1033 return False 1034 1035 # Separate the file name and path. 1036 path, file = os.path.split(file_path) 1037 1038 # Convert the structure reading args into lists. 1039 if read_mol and not isinstance(read_mol, list): 1040 read_mol = [read_mol] 1041 if set_mol_name and not isinstance(set_mol_name, list): 1042 set_mol_name = [set_mol_name] 1043 if read_model and not isinstance(read_model, list): 1044 read_model = [read_model] 1045 if set_model_num and not isinstance(set_model_num, list): 1046 set_model_num = [set_model_num] 1047 1048 # Loop over all models in the PDB file. 1049 model_index = 0 1050 orig_model_num = [] 1051 mol_conts = [] 1052 for model_num, model_records in self._parse_models_pdb(file_path): 1053 # Only load the desired model. 1054 if read_model and model_num not in read_model: 1055 continue 1056 1057 # Store the original model number. 1058 orig_model_num.append(model_num) 1059 1060 # Loop over the molecules of the model. 1061 mol_conts.append([]) 1062 mol_index = 0 1063 orig_mol_num = [] 1064 new_mol_name = [] 1065 for mol_num, mol_records in self._parse_mols(model_records): 1066 # Only load the desired model. 1067 if read_mol and mol_num not in read_mol: 1068 continue 1069 1070 # Set the target molecule name. 1071 if set_mol_name: 1072 new_mol_name.append(set_mol_name[mol_index]) 1073 else: 1074 # Number of structures already present for the model. 1075 num_struct = 0 1076 for model in self.structural_data: 1077 if not set_model_num or (model_index <= len(set_model_num) and set_model_num[model_index] == model.num): 1078 num_struct = len(model.mol) 1079 1080 # Set the name to the file name plus the structure number. 1081 new_mol_name.append(file_root(file) + '_mol' + repr(mol_num+num_struct)) 1082 1083 # Store the original mol number. 1084 orig_mol_num.append(mol_num) 1085 1086 # Generate the molecule container. 1087 mol = MolContainer() 1088 1089 # Fill the molecular data object. 1090 mol.fill_object_from_pdb(mol_records) 1091 1092 # Store the molecule container. 1093 mol_conts[model_index].append(mol) 1094 1095 # Increment the molecule index. 1096 mol_index = mol_index + 1 1097 1098 # Increment the model index. 1099 model_index = model_index + 1 1100 1101 # No data, so throw a warning and exit. 1102 if not len(mol_conts): 1103 warn(RelaxWarning("No structural data could be read from the file '%s'." % file_path)) 1104 return False 1105 1106 # Create the structural data data structures. 1107 self.pack_structs(mol_conts, orig_model_num=orig_model_num, set_model_num=set_model_num, orig_mol_num=orig_mol_num, set_mol_name=new_mol_name, file_name=file, file_path=path) 1108 1109 # Loading worked. 1110 return True
1111 1112
1113 - def load_xyz(self, file_path, read_mol=None, set_mol_name=None, read_model=None, set_model_num=None, verbosity=False):
1114 """Method for loading structures from a XYZ file. 1115 1116 @param file_path: The full path of the XYZ file. 1117 @type file_path: str 1118 @keyword read_mol: The molecule(s) to read from the file, independent of model. The 1119 molecules are determined differently by the different parsers, but 1120 are numbered consecutively from 1. If set to None, then all 1121 molecules will be loaded. 1122 @type read_mol: None, int, or list of int 1123 @keyword set_mol_name: Set the names of the molecules which are loaded. If set to None, 1124 then the molecules will be automatically labelled based on the file 1125 name or other information. 1126 @type set_mol_name: None, str, or list of str 1127 @keyword read_model: The XYZ model to extract from the file. If set to None, then all 1128 models will be loaded. 1129 @type read_model: None, int, or list of int 1130 @keyword set_model_num: Set the model number of the loaded molecule. If set to None, then 1131 the XYZ model numbers will be preserved, if they exist. 1132 @type set_model_num: None, int, or list of int 1133 @keyword verbosity: A flag which if True will cause messages to be printed. 1134 @type verbosity: bool 1135 @return: The status of the loading of the XYZ file. 1136 @rtype: bool 1137 """ 1138 1139 # Initial printout. 1140 if verbosity: 1141 print("\nInternal relax XYZ parser.") 1142 1143 # Test if the file exists. 1144 if not access(file_path, F_OK): 1145 # Exit indicating failure. 1146 return False 1147 1148 # Separate the file name and path. 1149 path, file = os.path.split(file_path) 1150 1151 # Convert the structure reading args into lists. 1152 if read_mol and not isinstance(read_mol, list): 1153 read_mol = [read_mol] 1154 if set_mol_name and not isinstance(set_mol_name, list): 1155 set_mol_name = [set_mol_name] 1156 if read_model and not isinstance(read_model, list): 1157 read_model = [read_model] 1158 if set_model_num and not isinstance(set_model_num, list): 1159 set_model_num = [set_model_num] 1160 1161 # Loop over all models in the XYZ file. 1162 mol_index=0 1163 model_index = 0 1164 xyz_model_increment = 0 1165 orig_model_num = [] 1166 mol_conts = [] 1167 orig_mol_num = [] 1168 new_mol_name = [] 1169 for model_records in self._parse_models_xyz(file_path): 1170 # Increment the xyz_model_increment 1171 xyz_model_increment = xyz_model_increment +1 1172 1173 # Only load the desired model. 1174 if read_model and xyz_model_increment not in read_model: 1175 continue 1176 1177 # Store the original model number. 1178 orig_model_num.append(model_index) 1179 1180 # Loop over the molecules of the model. 1181 if read_mol and mol_index not in read_mol: 1182 continue 1183 1184 # Set the target molecule name. 1185 if set_mol_name: 1186 new_mol_name.append(set_mol_name[mol_index]) 1187 else: 1188 if mol_index==0: 1189 #Set the name to the file name plus the structure number. 1190 new_mol_name.append(file_root(file) + '_mol' + repr(mol_index+1)) 1191 1192 # Store the original mol number. 1193 orig_mol_num.append(mol_index) 1194 1195 # Generate the molecule container. 1196 mol = MolContainer() 1197 1198 # Fill the molecular data object. 1199 mol.fill_object_from_xyz(model_records) 1200 1201 # Store the molecule container. 1202 mol_conts.append([]) 1203 mol_conts[model_index].append(mol) 1204 1205 # Increment the molecule index. 1206 mol_index = mol_index + 1 1207 1208 # Increment the model index. 1209 model_index = model_index + 1 1210 1211 orig_mol_num=[0] 1212 # Create the structural data data structures. 1213 self.pack_structs(mol_conts, orig_model_num=orig_model_num, set_model_num=set_model_num, orig_mol_num=orig_mol_num, set_mol_name=new_mol_name, file_name=file, file_path=path) 1214 1215 # Loading worked. 1216 return True
1217 1218
1219 - def rotate(self, R=None, origin=None, model=None, atom_id=None):
1220 """Rotate the structural information about the given origin. 1221 1222 @keyword R: The forwards rotation matrix. 1223 @type R: numpy 3D, rank-2 array 1224 @keyword origin: The origin of the rotation. 1225 @type origin: numpy 3D, rank-1 array 1226 @keyword model: The model to rotate. If None, all models will be rotated. 1227 @type model: int 1228 @keyword atom_id: The molecule, residue, and atom identifier string. Only atoms matching this selection will be used. 1229 @type atom_id: str or None 1230 """ 1231 1232 # Generate the selection object. 1233 sel_obj = None 1234 if atom_id: 1235 sel_obj = Selection(atom_id) 1236 1237 # Loop over the models. 1238 for model_cont in self.model_loop(model): 1239 # Loop over the molecules. 1240 for mol in model_cont.mol: 1241 # Skip non-matching molecules. 1242 if sel_obj and not sel_obj.contains_mol(mol.mol_name): 1243 continue 1244 1245 # Loop over the atoms. 1246 for i in range(len(mol.atom_num)): 1247 # Skip non-matching atoms. 1248 if sel_obj and not sel_obj.contains_spin(mol.atom_num[i], mol.atom_name[i], mol.res_num[i], mol.res_name[i], mol.mol_name): 1249 continue 1250 1251 # The origin to atom vector. 1252 vect = array([mol.x[i], mol.y[i], mol.z[i]], float64) - origin 1253 1254 # Rotation. 1255 rot_vect = dot(R, vect) 1256 1257 # The new position. 1258 pos = rot_vect + origin 1259 mol.x[i] = pos[0] 1260 mol.y[i] = pos[1] 1261 mol.z[i] = pos[2]
1262 1263
1264 - def translate(self, T=None, model=None, atom_id=None):
1265 """Displace the structural information by the given translation vector. 1266 1267 @keyword T: The translation vector. 1268 @type T: numpy 3D, rank-1 array 1269 @keyword model: The model to rotate. If None, all models will be rotated. 1270 @type model: int 1271 @keyword atom_id: The molecule, residue, and atom identifier string. Only atoms matching this selection will be used. 1272 @type atom_id: str or None 1273 """ 1274 1275 # Generate the selection object. 1276 sel_obj = None 1277 if atom_id: 1278 sel_obj = Selection(atom_id) 1279 1280 # Loop over the models. 1281 for model_cont in self.model_loop(model): 1282 # Loop over the molecules. 1283 for mol in model_cont.mol: 1284 # Skip non-matching molecules. 1285 if sel_obj and not sel_obj.contains_mol(mol.mol_name): 1286 continue 1287 1288 # Loop over the atoms. 1289 for i in range(len(mol.atom_num)): 1290 # Skip non-matching atoms. 1291 if sel_obj and not sel_obj.contains_spin(mol.atom_num[i], mol.atom_name[i], mol.res_num[i], mol.res_name[i], mol.mol_name): 1292 continue 1293 1294 # Translate. 1295 mol.x[i] = mol.x[i] + T[0] 1296 mol.y[i] = mol.y[i] + T[1] 1297 mol.z[i] = mol.z[i] + T[2]
1298 1299
1300 - def validate_models(self):
1301 """Check that the models are consistent with each other. 1302 1303 This checks that the primary structure is identical between the models. 1304 """ 1305 1306 # Print out. 1307 print("Validating models:") 1308 1309 # Loop over the models. 1310 for i in range(len(self.structural_data)): 1311 # Check the molecules. 1312 if len(self.structural_data[0].mol) != len(self.structural_data[i].mol): 1313 raise RelaxError("The number of molecules, %i, in model %i does not match the %i molecules of the first model." % (len(self.structural_data[i].mol), self.structural_data[i].num, len(self.structural_data[0].mol))) 1314 1315 # Loop over the molecules. 1316 for j in range(len(self.structural_data[i].mol)): 1317 # Alias the molecules. 1318 mol = self.structural_data[i].mol[j] 1319 mol_ref = self.structural_data[0].mol[j] 1320 1321 # Check the names. 1322 if mol.mol_name != mol_ref.mol_name: 1323 raise RelaxError("The molecule name '%s' of model %i does not match the name '%s' of the first model." % (mol.mol_name, self.structural_data[i].num, mol_ref.mol_name)) 1324 1325 # Loop over the atoms. 1326 for k in range(len(mol.atom_name)): 1327 # Create pseudo-pdb formatted records (with no atomic coordinates). 1328 atom = "%-6s%5s %4s%1s%3s %1s%4s%1s %8s%8s%8s%6.2f%6.2f %4s%2s%2s" % ('ATOM', mol.atom_num[k], self._translate(mol.atom_name[k]), '', self._translate(mol.res_name[k]), self._translate(mol.chain_id[k]), self._translate(mol.res_num[k]), '', '#', '#', '#', 1.0, 0, self._translate(mol.seg_id[k]), self._translate(mol.element[k]), '') 1329 atom_ref = "%-6s%5s %4s%1s%3s %1s%4s%1s %8s%8s%8s%6.2f%6.2f %4s%2s%2s" % ('ATOM', mol_ref.atom_num[k], self._translate(mol_ref.atom_name[k]), '', self._translate(mol_ref.res_name[k]), self._translate(mol_ref.chain_id[k]), self._translate(mol_ref.res_num[k]), '', '#', '#', '#', 1.0, 0, self._translate(mol_ref.seg_id[k]), self._translate(mol_ref.element[k]), '') 1330 1331 # Check the atom info. 1332 if atom != atom_ref: 1333 print(atom) 1334 print(atom_ref) 1335 raise RelaxError("The atoms of model %i do not match the first model." % self.structural_data[i].num) 1336 1337 # Final printout. 1338 print("\tAll models are consistent")
1339 1340
1341 - def write_pdb(self, file, model_num=None):
1342 """Method for the creation of a PDB file from the structural data. 1343 1344 A number of PDB records including HET, HETNAM, FORMUL, HETATM, TER, CONECT, MASTER, and END 1345 are created. To create the non-standard residue records HET, HETNAM, and FORMUL, the data 1346 structure 'het_data' is created. It is an array of arrays where the first dimension 1347 corresponds to a different residue and the second dimension has the elements: 1348 1349 0. Residue number. 1350 1. Residue name. 1351 2. Chain ID. 1352 3. Total number of atoms in the residue. 1353 4. Number of H atoms in the residue. 1354 5. Number of C atoms in the residue. 1355 1356 1357 @param file: The PDB file object. This object must be writable. 1358 @type file: file object 1359 @keyword model_num: The model to place into the PDB file. If not supplied, then all 1360 models will be placed into the file. 1361 @type model_num: None or int 1362 """ 1363 1364 # Validate the structural data. 1365 self.validate() 1366 1367 # Initialise record counts. 1368 num_hetatm = 0 1369 num_atom = 0 1370 num_ter = 0 1371 num_conect = 0 1372 1373 # Print out. 1374 print("\nCreating the PDB records\n") 1375 1376 # Write some initial remarks. 1377 print("REMARK") 1378 file.write("REMARK 4 THIS FILE COMPLIES WITH FORMAT V. 3.1, 1-AUG-2007\n") 1379 file.write("REMARK 40 CREATED BY RELAX (HTTP://NMR-RELAX.COM)\n") 1380 num_remark = 2 1381 1382 # Determine if model records will be created. 1383 model_records = False 1384 for model in self.model_loop(): 1385 if hasattr(model, 'num') and model.num != None: 1386 model_records = True 1387 1388 1389 #################### 1390 # Hetrogen section # 1391 #################### 1392 1393 # Initialise the hetrogen info array. 1394 het_data = [] 1395 het_data_coll = [] 1396 1397 # Loop over the molecules of the first model. 1398 index = 0 1399 for mol in self.structural_data[0].mol: 1400 # Check the validity of the data. 1401 self._validate_data_arrays(mol) 1402 1403 # Append an empty array for this molecule. 1404 het_data.append([]) 1405 1406 # Collect the non-standard residue info. 1407 for i in range(len(mol.atom_name)): 1408 # Skip non-HETATM records and HETATM records with no residue info. 1409 if mol.pdb_record[i] != 'HETATM' or mol.res_name[i] == None: 1410 continue 1411 1412 # If the residue is not already stored initialise a new het_data element. 1413 # (residue number, residue name, chain ID, number of atoms, atom count array). 1414 if not het_data[index] or not mol.res_num[i] == het_data[index][-1][0]: 1415 het_data[index].append([mol.res_num[i], mol.res_name[i], mol.chain_id[i], 0, []]) 1416 1417 # Catch missing chain_ids. 1418 if het_data[index][-1][2] == None: 1419 het_data[index][-1][2] = '' 1420 1421 # Total atom count. 1422 het_data[index][-1][3] = het_data[index][-1][3] + 1 1423 1424 # Find if the atom has already a count entry. 1425 entry = False 1426 for j in range(len(het_data[index][-1][4])): 1427 if mol.element[i] == het_data[index][-1][4][j][0]: 1428 entry = True 1429 1430 # Create a new specific atom count entry. 1431 if not entry: 1432 het_data[index][-1][4].append([mol.element[i], 0]) 1433 1434 # Increment the specific atom count. 1435 for j in range(len(het_data[index][-1][4])): 1436 if mol.element[i] == het_data[index][-1][4][j][0]: 1437 het_data[index][-1][4][j][1] = het_data[index][-1][4][j][1] + 1 1438 1439 # Create the collective hetrogen info data structure. 1440 for i in range(len(het_data[index])): 1441 # Find the entry in the collective structure. 1442 found = False 1443 for j in range(len(het_data_coll)): 1444 # Matching residue numbers. 1445 if het_data[index][i][0] == het_data_coll[j][0]: 1446 # Change the flag. 1447 found = True 1448 1449 # The checks. 1450 if het_data_coll[j][1] != het_data[index][i][1]: 1451 raise RelaxError("The " + repr(het_data[index][i][1]) + " residue name of hetrogen " + repr(het_data[index][i][0]) + " " + het_data[index][i][1] + " of structure " + repr(index) + " does not match the " + repr(het_data_coll[j][1]) + " name of the previous structures.") 1452 1453 elif het_data_coll[j][2] != het_data[index][i][2]: 1454 raise RelaxError("The hetrogen chain id " + repr(het_data[index][i][2]) + " does not match " + repr(het_data_coll[j][2]) + " of residue " + repr(het_data_coll[j][0]) + " " + het_data_coll[j][1] + " of the previous structures.") 1455 1456 elif het_data_coll[j][3] != het_data[index][i][3]: 1457 raise RelaxError("The " + repr(het_data[index][i][3]) + " atoms of hetrogen " + repr(het_data_coll[j][0]) + " " + het_data_coll[j][1] + " of structure " + repr(index) + " does not match the " + repr(het_data_coll[j][3]) + " of the previous structures.") 1458 1459 elif het_data_coll[j][4] != het_data[index][i][4]: 1460 raise RelaxError("The atom counts " + repr(het_data[index][i][4]) + " for the hetrogen residue " + repr(het_data_coll[j][0]) + " " + het_data_coll[j][1] + " of structure " + repr(index) + " do not match the counts " + repr(het_data_coll[j][4]) + " of the previous structures.") 1461 1462 # If there is no match, add the new residue to the collective. 1463 if not found: 1464 het_data_coll.append(het_data[index][i]) 1465 1466 # Increment the molecule index. 1467 index = index + 1 1468 1469 1470 # The HET records. 1471 ################## 1472 1473 # Print out. 1474 print("HET") 1475 1476 # Write the HET records. 1477 for het in het_data_coll: 1478 file.write("%-6s %3s %1s%4s%1s %5s %-40s\n" % ('HET', het[2], het[1], het[0], '', het[3], '')) 1479 1480 1481 # The HETNAM records. 1482 ##################### 1483 1484 # Print out. 1485 print("HETNAM") 1486 1487 # Loop over the non-standard residues. 1488 residues = [] 1489 for het in het_data_coll: 1490 # Test if the residue HETNAM record as already been written (otherwise store its name). 1491 if het[1] in residues: 1492 continue 1493 else: 1494 residues.append(het[1]) 1495 1496 # Get the chemical name. 1497 chemical_name = self._get_chemical_name(het[1]) 1498 if not chemical_name: 1499 chemical_name = 'Unknown' 1500 1501 # Write the HETNAM records. 1502 file.write("%-6s %2s %3s %-55s\n" % ('HETNAM', '', het[1], chemical_name)) 1503 1504 1505 # The FORMUL records. 1506 ##################### 1507 1508 # Print out. 1509 print("FORMUL") 1510 1511 # Loop over the non-standard residues and generate and write the chemical formula. 1512 residues = [] 1513 for het in het_data_coll: 1514 # Test if the residue HETNAM record as already been written (otherwise store its name). 1515 if het[1] in residues: 1516 continue 1517 else: 1518 residues.append(het[1]) 1519 1520 # Initialise the chemical formula. 1521 formula = '' 1522 1523 # Loop over the atoms. 1524 for atom_count in het[4]: 1525 formula = formula + atom_count[0] + repr(atom_count[1]) 1526 1527 # The FORMUL record (chemical formula). 1528 file.write("%-6s %2s %3s %2s%1s%-51s\n" % ('FORMUL', het[0], het[1], '', '', formula)) 1529 1530 1531 ###################### 1532 # Coordinate section # 1533 ###################### 1534 1535 # Loop over the models. 1536 for model in self.model_loop(model_num): 1537 # MODEL record, for multiple models. 1538 #################################### 1539 1540 if model_records: 1541 # Print out. 1542 print("\nMODEL %s" % model.num) 1543 1544 # Write the model record. 1545 file.write("%-6s %4i\n" % ('MODEL', model.num)) 1546 1547 1548 # Add the atomic coordinate records (ATOM, HETATM, and TER). 1549 ############################################################ 1550 1551 # Loop over the molecules. 1552 for mol in model.mol: 1553 # Print out. 1554 print("ATOM, HETATM, TER") 1555 1556 # Loop over the atomic data. 1557 atom_record = False 1558 for i in range(len(mol.atom_name)): 1559 # Write the ATOM record. 1560 if mol.pdb_record[i] in [None, 'ATOM']: 1561 atom_record = True 1562 1563 # The atom number, if missing. 1564 atom_num = mol.atom_num[i] 1565 if atom_num == None: 1566 atom_num = i + 1 1567 1568 # Write out. 1569 file.write("%-6s%5s %4s%1s%3s %1s%4s%1s %8.3f%8.3f%8.3f%6.2f%6.2f %4s%2s%2s\n" % ('ATOM', atom_num, self._translate(mol.atom_name[i]), '', self._translate(mol.res_name[i]), self._translate(mol.chain_id[i]), self._translate(mol.res_num[i]), '', self._translate(mol.x[i], 'float'), self._translate(mol.y[i], 'float'), self._translate(mol.z[i], 'float'), 1.0, 0, self._translate(mol.seg_id[i]), self._translate(mol.element[i]), '')) 1570 num_atom = num_atom + 1 1571 1572 # Info for the TER record. 1573 ter_num = atom_num + 1 1574 ter_name = mol.res_name[i] 1575 ter_chain_id = mol.chain_id[i] 1576 ter_res_num = mol.res_num[i] 1577 1578 # Finish the ATOM section with the TER record. 1579 if atom_record: 1580 file.write("%-6s%5s %3s %1s%4s%1s\n" % ('TER', ter_num, self._translate(ter_name), self._translate(ter_chain_id), self._translate(ter_res_num), '')) 1581 num_ter = num_ter + 1 1582 1583 # Loop over the atomic data. 1584 count_shift = False 1585 for i in range(len(mol.atom_name)): 1586 # Write the HETATM record. 1587 if mol.pdb_record[i] == 'HETATM': 1588 # The atom number, if missing. 1589 atom_num = mol.atom_num[i] 1590 if atom_num == None: 1591 atom_num = i + 1 1592 1593 # Increment the atom number if a TER record was created. 1594 if atom_record and atom_num == ter_num: 1595 count_shift = True 1596 if atom_record and count_shift: 1597 atom_num += 1 1598 1599 # Write out. 1600 file.write("%-6s%5s %4s%1s%3s %1s%4s%1s %8.3f%8.3f%8.3f%6.2f%6.2f %4s%2s%2s\n" % ('HETATM', atom_num, self._translate(mol.atom_name[i]), '', self._translate(mol.res_name[i]), self._translate(mol.chain_id[i]), self._translate(mol.res_num[i]), '', self._translate(mol.x[i], 'float'), self._translate(mol.y[i], 'float'), self._translate(mol.z[i], 'float'), 1.0, 0, self._translate(mol.seg_id[i]), self._translate(mol.element[i]), '')) 1601 num_hetatm = num_hetatm + 1 1602 1603 1604 # ENDMDL record, for multiple structures. 1605 ######################################## 1606 1607 if model_records: 1608 # Print out. 1609 print("ENDMDL") 1610 1611 # Write the model record. 1612 file.write("%-6s\n" % 'ENDMDL') 1613 1614 1615 # Create the CONECT records. 1616 ############################ 1617 1618 # Print out. 1619 print("CONECT") 1620 1621 # Loop over the molecules of the first model. 1622 for mol in self.structural_data[0].mol: 1623 # Loop over the atoms. 1624 for i in range(len(mol.atom_name)): 1625 # No bonded atoms, hence no CONECT record is required. 1626 if not len(mol.bonded[i]): 1627 continue 1628 1629 # Initialise some data structures. 1630 flush = 0 1631 bonded_index = 0 1632 bonded = ['', '', '', ''] 1633 1634 # Loop over the bonded atoms. 1635 for j in range(len(mol.bonded[i])): 1636 # End of the array, hence create the CONECT record in this iteration. 1637 if j == len(mol.bonded[i])-1: 1638 flush = True 1639 1640 # Only four covalently bonded atoms allowed in one CONECT record. 1641 if bonded_index == 3: 1642 flush = True 1643 1644 # Get the bonded atom index. 1645 bonded[bonded_index] = mol.bonded[i][j] 1646 1647 # Increment the bonded_index value. 1648 bonded_index = bonded_index + 1 1649 1650 # Generate the CONECT record and increment the counter. 1651 if flush: 1652 # Convert the atom indices to atom numbers. 1653 for k in range(4): 1654 if bonded[k] != '': 1655 if mol.atom_num[bonded[k]] != None: 1656 bonded[k] = mol.atom_num[bonded[k]] 1657 else: 1658 bonded[k] = bonded[k] + 1 1659 1660 # Write the CONECT record. 1661 file.write("%-6s%5s%5s%5s%5s%5s%5s%5s%5s%5s%5s%5s\n" % ('CONECT', i+1, bonded[0], bonded[1], bonded[2], bonded[3], '', '', '', '', '', '')) 1662 1663 # Reset the flush flag, the bonded atom count, and the bonded atom names. 1664 flush = False 1665 bonded_index = 0 1666 bonded = ['', '', '', ''] 1667 1668 # Increment the CONECT record count. 1669 num_conect = num_conect + 1 1670 1671 1672 1673 # MASTER record. 1674 ################ 1675 1676 # Print out. 1677 print("\nMASTER") 1678 1679 # Write the MASTER record. 1680 file.write("%-6s %5s%5s%5s%5s%5s%5s%5s%5s%5s%5s%5s%5s\n" % ('MASTER', 0, 0, len(het_data_coll), 0, 0, 0, 0, 0, num_atom+num_hetatm, num_ter, num_conect, 0)) 1681 1682 1683 # END. 1684 ###### 1685 1686 # Print out. 1687 print("END") 1688 1689 # Write the END record. 1690 file.write("END\n")
1691 1692
1693 -class MolContainer:
1694 """The container for the molecular information. 1695 1696 The structural data object for this class is a container possessing a number of different arrays 1697 corresponding to different structural information. These objects include: 1698 1699 - atom_num: The atom name. 1700 - atom_name: The atom name. 1701 - bonded: Each element an array of bonded atom indices. 1702 - chain_id: The chain ID. 1703 - element: The element symbol. 1704 - pdb_record: The optional PDB record name (one of ATOM, HETATM, or TER). 1705 - res_name: The residue name. 1706 - res_num: The residue number. 1707 - seg_id: The segment ID. 1708 - x: The x coordinate of the atom. 1709 - y: The y coordinate of the atom. 1710 - z: The z coordinate of the atom. 1711 1712 All arrays should be of equal length so that an atom index can retrieve all the corresponding 1713 data. Only the atom identification string is compulsory, all other arrays can contain None. 1714 """ 1715 1716
1717 - def __init__(self):
1718 """Initialise the molecular container.""" 1719 1720 # The atom num (array of int). 1721 self.atom_num = [] 1722 1723 # The atom name (array of str). 1724 self.atom_name = [] 1725 1726 # The bonded atom indices (array of arrays of int). 1727 self.bonded = [] 1728 1729 # The chain ID (array of str). 1730 self.chain_id = [] 1731 1732 # The element symbol (array of str). 1733 self.element = [] 1734 1735 # The optional PDB record name (array of str). 1736 self.pdb_record = [] 1737 1738 # The residue name (array of str). 1739 self.res_name = [] 1740 1741 # The residue number (array of int). 1742 self.res_num = [] 1743 1744 # The segment ID (array of int). 1745 self.seg_id = [] 1746 1747 # The x coordinate (array of float). 1748 self.x = [] 1749 1750 # The y coordinate (array of float). 1751 self.y = [] 1752 1753 # The z coordinate (array of float). 1754 self.z = []
1755 1756
1757 - def _atom_index(self, atom_num):
1758 """Find the atom index corresponding to the given atom number. 1759 1760 @param atom_num: The atom number to find the index of. 1761 @type atom_num: int 1762 @return: The atom index corresponding to the atom. 1763 @rtype: int 1764 """ 1765 1766 # Loop over the atoms. 1767 for j in range(len(self.atom_num)): 1768 # Return the index. 1769 if self.atom_num[j] == atom_num: 1770 return j 1771 1772 # Should not be here, the PDB connect records are incorrect. 1773 warn(RelaxWarning("The atom number " + repr(atom_num) + " from the CONECT record cannot be found within the ATOM and HETATM records."))
1774 1775
1776 - def _det_pdb_element(self, atom_name):
1777 """Try to determine the element from the PDB atom name. 1778 1779 @param atom_name: The PDB atom name. 1780 @type atom_name: str 1781 @return: The element name, or None if unsuccessful. 1782 @rtype: str or None 1783 """ 1784 1785 # Strip away the "'" character (for RNA, etc.). 1786 element = atom_name.strip("'") 1787 1788 # Strip away atom numbering, from the front and end. 1789 element = element.strip(digits) 1790 1791 # Amino acid atom translation table (note, numbers have been stripped already!). 1792 table = {'C': ['CA', 'CB', 'CG', 'CD', 'CE', 'CH', 'CZ'], 1793 'N': ['ND', 'NE', 'NH', 'NZ'], 1794 'H': ['HA', 'HB', 'HG', 'HD', 'HE', 'HH', 'HT', 'HZ'], 1795 'O': ['OG', 'OD', 'OE', 'OH', 'OT'], 1796 'S': ['SD', 'SG'] 1797 } 1798 1799 # Translate amino acids. 1800 for key in list(table.keys()): 1801 if element in table[key]: 1802 element = key 1803 break 1804 1805 # Allowed element list. 1806 elements = ['H', 'C', 'N', 'O', 'F', 'P', 'S'] 1807 1808 # Return the element, if in the list. 1809 if element in elements: 1810 return element 1811 1812 # Else, throw a warning. 1813 warn(RelaxWarning("Cannot determine the element associated with atom '%s'." % atom_name))
1814 1815
1816 - def _parse_pdb_record(self, record):
1817 """Parse the PDB record string and return an array of the corresponding atomic information. 1818 1819 The format of the ATOM and HETATM records is:: 1820 __________________________________________________________________________________________ 1821 | | | | | 1822 | Columns | Data type | Field | Definition | 1823 |_________|______________|______________|________________________________________________| 1824 | | | | | 1825 | 1 - 6 | Record name | "ATOM" | | 1826 | 7 - 11 | Integer | serial | Atom serial number. | 1827 | 13 - 16 | Atom | name | Atom name. | 1828 | 17 | Character | altLoc | Alternate location indicator. | 1829 | 18 - 20 | Residue name | resName | Residue name. | 1830 | 22 | Character | chainID | Chain identifier. | 1831 | 23 - 26 | Integer | resSeq | Residue sequence number. | 1832 | 27 | AChar | iCode | Code for insertion of residues. | 1833 | 31 - 38 | Real(8.3) | x | Orthogonal coordinates for X in Angstroms. | 1834 | 39 - 46 | Real(8.3) | y | Orthogonal coordinates for Y in Angstroms. | 1835 | 47 - 54 | Real(8.3) | z | Orthogonal coordinates for Z in Angstroms. | 1836 | 55 - 60 | Real(6.2) | occupancy | Occupancy. | 1837 | 61 - 66 | Real(6.2) | tempFactor | Temperature factor. | 1838 | 73 - 76 | LString(4) | segID | Segment identifier, left-justified. | 1839 | 77 - 78 | LString(2) | element | Element symbol, right-justified. | 1840 | 79 - 80 | LString(2) | charge | Charge on the atom. | 1841 |_________|______________|______________|________________________________________________| 1842 1843 1844 The format of the TER record is:: 1845 __________________________________________________________________________________________ 1846 | | | | | 1847 | Columns | Data type | Field | Definition | 1848 |_________|______________|______________|________________________________________________| 1849 | | | | | 1850 | 1 - 6 | Record name | "TER " | | 1851 | 7 - 11 | Integer | serial | Serial number. | 1852 | 18 - 20 | Residue name | resName | Residue name. | 1853 | 22 | Character | chainID | Chain identifier. | 1854 | 23 - 26 | Integer | resSeq | Residue sequence number. | 1855 | 27 | AChar | iCode | Insertion code. | 1856 |_________|______________|______________|________________________________________________| 1857 1858 1859 The format of the CONECT record is:: 1860 __________________________________________________________________________________________ 1861 | | | | | 1862 | Columns | Data type | Field | Definition | 1863 |_________|______________|______________|________________________________________________| 1864 | | | | | 1865 | 1 - 6 | Record name | "CONECT" | | 1866 | 7 - 11 | Integer | serial | Atom serial number | 1867 | 12 - 16 | Integer | serial | Serial number of bonded atom | 1868 | 17 - 21 | Integer | serial | Serial number of bonded atom | 1869 | 22 - 26 | Integer | serial | Serial number of bonded atom | 1870 | 27 - 31 | Integer | serial | Serial number of bonded atom | 1871 | 32 - 36 | Integer | serial | Serial number of hydrogen bonded atom | 1872 | 37 - 41 | Integer | serial | Serial number of hydrogen bonded atom | 1873 | 42 - 46 | Integer | serial | Serial number of salt bridged atom | 1874 | 47 - 51 | Integer | serial | Serial number of hydrogen bonded atom | 1875 | 52 - 56 | Integer | serial | Serial number of hydrogen bonded atom | 1876 | 57 - 61 | Integer | serial | Serial number of salt bridged atom | 1877 |_________|______________|______________|________________________________________________| 1878 1879 1880 @param record: The single line PDB record. 1881 @type record: str 1882 @return: The list of atomic information, each element corresponding to the PDB fields 1883 as defined in "Protein Data Bank Contents Guide: Atomic Coordinate Entry 1884 Format Description" version 2.1 (draft), October 25, 1996. 1885 @rtype: list of str 1886 """ 1887 1888 # Initialise. 1889 fields = [] 1890 1891 # ATOM and HETATM records. 1892 if search('^ATOM', record) or search('^HETATM', record): 1893 # Split up the record. 1894 fields.append(record[0:6]) 1895 fields.append(record[6:11]) 1896 fields.append(record[12:16]) 1897 fields.append(record[16]) 1898 fields.append(record[17:20]) 1899 fields.append(record[21]) 1900 fields.append(record[22:26]) 1901 fields.append(record[26]) 1902 fields.append(record[30:38]) 1903 fields.append(record[38:46]) 1904 fields.append(record[46:54]) 1905 fields.append(record[54:60]) 1906 fields.append(record[60:66]) 1907 fields.append(record[72:76]) 1908 fields.append(record[76:78]) 1909 fields.append(record[78:80]) 1910 1911 # Loop over the fields. 1912 for i in range(len(fields)): 1913 # Strip all whitespace. 1914 fields[i] = fields[i].strip() 1915 1916 # Replace nothingness with None. 1917 if fields[i] == '': 1918 fields[i] = None 1919 1920 # Convert strings to numbers. 1921 if fields[1]: 1922 fields[1] = int(fields[1]) 1923 if fields[6]: 1924 fields[6] = int(fields[6]) 1925 if fields[8]: 1926 fields[8] = float(fields[8]) 1927 if fields[9]: 1928 fields[9] = float(fields[9]) 1929 if fields[10]: 1930 fields[10] = float(fields[10]) 1931 if fields[11]: 1932 fields[11] = float(fields[11]) 1933 if fields[12]: 1934 fields[12] = float(fields[12]) 1935 1936 # CONECT records. 1937 if search('^CONECT', record): 1938 # Split up the record. 1939 fields.append(record[0:6]) 1940 fields.append(record[6:11]) 1941 fields.append(record[11:16]) 1942 fields.append(record[16:21]) 1943 fields.append(record[21:26]) 1944 fields.append(record[26:31]) 1945 1946 # Loop over the fields. 1947 for i in range(len(fields)): 1948 # Strip all whitespace. 1949 fields[i] = fields[i].strip() 1950 1951 # Replace nothingness with None. 1952 if fields[i] == '': 1953 fields[i] = None 1954 1955 # Convert strings to numbers. 1956 if fields[1]: 1957 fields[1] = int(fields[1]) 1958 if fields[2]: 1959 fields[2] = int(fields[2]) 1960 if fields[3]: 1961 fields[3] = int(fields[3]) 1962 if fields[4]: 1963 fields[4] = int(fields[4]) 1964 if fields[5]: 1965 fields[5] = int(fields[5]) 1966 1967 # Return the atomic info. 1968 return fields
1969 1970
1971 - def _parse_xyz_record(self, record):
1972 """Parse the XYZ record string and return an array of the corresponding atomic information. 1973 1974 The format of the XYZ records is:: 1975 __________________________________________________________________________________________ 1976 | | | | | 1977 | Columns | Data type | Field | Definition | 1978 |_________|______________|______________|________________________________________________| 1979 | | | | | 1980 | 1 | String | element | | 1981 | 2 | Real | x | Orthogonal coordinates for X in Angstroms | 1982 | 3 | Real | y | Orthogonal coordinates for Y in Angstroms | 1983 | 4 | Real | z | Orthogonal coordinates for Z in Angstroms | 1984 |_________|______________|______________|________________________________________________| 1985 1986 1987 @param record: The single line PDB record. 1988 @type record: str 1989 @return: The list of atomic information 1990 @rtype: list of str 1991 """ 1992 1993 # Initialise. 1994 fields = [] 1995 word = record.split() 1996 1997 # ATOM and HETATM records. 1998 if len(word)==4: 1999 # Split up the record. 2000 fields.append(word[0]) 2001 fields.append(word[1]) 2002 fields.append(word[2]) 2003 fields.append(word[3]) 2004 2005 # Loop over the fields. 2006 for i in range(len(fields)): 2007 # Strip all whitespace. 2008 fields[i] = fields[i].strip() 2009 2010 # Replace nothingness with None. 2011 if fields[i] == '': 2012 fields[i] = None 2013 2014 # Convert strings to numbers. 2015 if fields[1]: 2016 fields[1] = float(fields[1]) 2017 if fields[2]: 2018 fields[2] = float(fields[2]) 2019 if fields[3]: 2020 fields[3] = float(fields[3]) 2021 2022 # Return the atomic info. 2023 return fields
2024 2025
2026 - def atom_add(self, 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):
2027 """Method for adding an atom to the structural data object. 2028 2029 This method will create the key-value pair for the given atom. 2030 2031 2032 @keyword atom_name: The atom name, e.g. 'H1'. 2033 @type atom_name: str or None 2034 @keyword res_name: The residue name. 2035 @type res_name: str or None 2036 @keyword res_num: The residue number. 2037 @type res_num: int or None 2038 @keyword pos: The position vector of coordinates. 2039 @type pos: list (length = 3) 2040 @keyword element: The element symbol. 2041 @type element: str or None 2042 @keyword atom_num: The atom number. 2043 @type atom_num: int or None 2044 @keyword chain_id: The chain identifier. 2045 @type chain_id: str or None 2046 @keyword segment_id: The segment identifier. 2047 @type segment_id: str or None 2048 @keyword pdb_record: The optional PDB record name, e.g. 'ATOM' or 'HETATM'. 2049 @type pdb_record: str or None 2050 """ 2051 2052 # Append to all the arrays. 2053 self.atom_num.append(atom_num) 2054 self.atom_name.append(atom_name) 2055 self.bonded.append([]) 2056 self.chain_id.append(chain_id) 2057 self.element.append(element) 2058 self.pdb_record.append(pdb_record) 2059 self.res_name.append(res_name) 2060 self.res_num.append(res_num) 2061 self.seg_id.append(segment_id) 2062 self.x.append(pos[0]) 2063 self.y.append(pos[1]) 2064 self.z.append(pos[2])
2065 2066
2067 - def atom_connect(self, index1=None, index2=None):
2068 """Method for connecting two atoms within the data structure object. 2069 2070 This method will append index2 to the array at bonded[index1] and vice versa. 2071 2072 2073 @keyword index1: The index of the first atom. 2074 @type index1: int 2075 @keyword index2: The index of the second atom. 2076 @type index2: int 2077 """ 2078 2079 # Update the bonded array structure, if necessary. 2080 if index2 not in self.bonded[index1]: 2081 self.bonded[index1].append(index2) 2082 if index1 not in self.bonded[index2]: 2083 self.bonded[index2].append(index1)
2084 2085
2086 - def fill_object_from_pdb(self, records):
2087 """Method for generating a complete Structure_container object from the given PDB records. 2088 2089 @param records: A list of structural PDB records. 2090 @type records: list of str 2091 """ 2092 2093 # Loop over the records. 2094 for record in records: 2095 # Parse the record. 2096 record = self._parse_pdb_record(record) 2097 2098 # Nothing to do. 2099 if not record: 2100 continue 2101 2102 # Add the atom. 2103 if record[0] == 'ATOM' or record[0] == 'HETATM': 2104 # Attempt at determining the element, if missing. 2105 element = record[14] 2106 if not element: 2107 element = self._det_pdb_element(record[2]) 2108 2109 # Add. 2110 self.atom_add(pdb_record=record[0], atom_num=record[1], atom_name=record[2], res_name=record[4], chain_id=record[5], res_num=record[6], pos=[record[8], record[9], record[10]], segment_id=record[13], element=element) 2111 2112 # Connect atoms. 2113 if record[0] == 'CONECT': 2114 # Loop over the atoms of the record. 2115 for i in range(len(record)-2): 2116 # Skip if there is no record. 2117 if record[i+2] == None: 2118 continue 2119 2120 # Make the connection. 2121 self.atom_connect(index1=self._atom_index(record[1]), index2=self._atom_index(record[i+2]))
2122 2123
2124 - def fill_object_from_xyz(self, records):
2125 """Method for generating a complete Structure_container object from the given xyz records. 2126 2127 @param records: A list of structural xyz records. 2128 @type records: list of str 2129 """ 2130 2131 # initialisation for atom number 2132 atom_number = 1 2133 2134 # Loop over the records. 2135 for record in records: 2136 # Parse the record. 2137 record = self._parse_xyz_record(record) 2138 2139 # Nothing to do. 2140 if not record: 2141 continue 2142 2143 # Add the atom. 2144 if len(record) == 4: 2145 # Attempt at determining the element, if missing. 2146 element = record[0] 2147 if not element: 2148 element = self._det_pdb_element(record[2]) 2149 2150 # Add. 2151 self.atom_add(atom_num=atom_number, pos=[record[1], record[2], record[3]], element=element) 2152 2153 # Increment of atom number 2154 atom_number = atom_number + 1
2155 2156
2157 - def from_xml(self, mol_node, file_version=1):
2158 """Recreate the MolContainer from the XML molecule node. 2159 2160 @param mol_node: The molecule XML node. 2161 @type mol_node: xml.dom.minicompat.NodeList instance 2162 @keyword file_version: The relax XML version of the XML file. 2163 @type file_version: int 2164 """ 2165 2166 # Recreate the current molecule container. 2167 xml_to_object(mol_node, self, file_version=file_version)
2168 2169
2170 - def is_empty(self):
2171 """Check if the container is empty.""" 2172 2173 # Set attributes. 2174 if hasattr(self, 'mol_name'): return False 2175 if hasattr(self, 'file_name'): return False 2176 if hasattr(self, 'file_path'): return False 2177 if hasattr(self, 'file_mol_num'): return False 2178 if hasattr(self, 'file_model'): return False 2179 2180 # Internal data structures. 2181 if not self.atom_num == []: return False 2182 if not self.atom_name == []: return False 2183 if not self.bonded == []: return False 2184 if not self.chain_id == []: return False 2185 if not self.element == []: return False 2186 if not self.pdb_record == []: return False 2187 if not self.res_name == []: return False 2188 if not self.res_num == []: return False 2189 if not self.seg_id == []: return False 2190 if not self.x == []: return False 2191 if not self.y == []: return False 2192 if not self.z == []: return False 2193 2194 # Ok, now this thing must be empty. 2195 return True
2196 2197
2198 - def to_xml(self, doc, element):
2199 """Create XML elements for the contents of this molecule container. 2200 2201 @param doc: The XML document object. 2202 @type doc: xml.dom.minidom.Document instance 2203 @param element: The element to add the molecule XML elements to. 2204 @type element: XML element object 2205 """ 2206 2207 # Create an XML element for this molecule and add it to the higher level element. 2208 mol_element = doc.createElement('mol_cont') 2209 element.appendChild(mol_element) 2210 2211 # Set the molecule attributes. 2212 mol_element.setAttribute('desc', 'Molecule container') 2213 mol_element.setAttribute('name', str(self.mol_name)) 2214 2215 # Add all simple python objects within the MolContainer to the XML element. 2216 fill_object_contents(doc, mol_element, object=self, blacklist=list(self.__class__.__dict__.keys()))
2217