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

Source Code for Module lib.structure.internal.object

   1  ############################################################################### 
   2  #                                                                             # 
   3  # Copyright (C) 2003-2014 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 copy import deepcopy 
  27  from numpy import array, dot, float64, linalg, zeros 
  28  import os 
  29  from os import F_OK, access, curdir, sep 
  30  from os.path import abspath 
  31  from re import search 
  32  from time import asctime 
  33  from warnings import warn 
  34   
  35  # relax module imports. 
  36  from lib import regex 
  37  from lib.check_types import is_float 
  38  from lib.errors import RelaxError, RelaxNoneIntError, RelaxNoPdbError 
  39  from lib.io import file_root, open_read_file 
  40  from lib.selection import Selection 
  41  from lib.structure import pdb_read, pdb_write 
  42  from lib.structure.internal.displacements import Displacements 
  43  from lib.structure.internal.models import ModelList 
  44  from lib.structure.internal.molecules import MolContainer 
  45  from lib.warnings import RelaxWarning 
  46  from lib.xml import object_to_xml, xml_to_object 
  47   
  48   
  49  # Module variables. 
  50  CHAIN_ID_LIST = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789abcdefghijklmnopqrstuvwxyz' 
  51  RELAX_VERSION = None 
  52   
  53   
54 -class Internal:
55 """The internal relax structural data object. 56 57 The structural data object for this class is a container possessing a number of different arrays corresponding to different structural information. These objects are described in the structural container docstring. 58 """ 59
60 - def __init__(self):
61 """Initialise the structural object.""" 62 63 # Initialise the empty model list. 64 self.structural_data = ModelList()
65 66
67 - def _bonded_atom(self, attached_atom, index, mol):
68 """Find the atom named attached_atom directly bonded to the atom located at the index. 69 70 @param attached_atom: The name of the attached atom to return. 71 @type attached_atom: str 72 @param index: The index of the atom which the attached atom is attached to. 73 @type index: int 74 @param mol: The molecule container. 75 @type mol: MolContainer instance 76 @return: A tuple of information about the bonded atom. 77 @rtype: tuple consisting of the atom number (int), atom name (str), element name (str), and atomic position (Numeric array of len 3) 78 """ 79 80 # Init. 81 bonded_found = False 82 83 # No bonded atoms, so determine the connectivities. 84 if not mol.bonded[index]: 85 # Determine the molecule type if needed. 86 if not hasattr(mol, 'type'): 87 self._mol_type(mol) 88 89 # Protein. 90 if mol.type == 'protein': 91 self._protein_connect(mol) 92 93 # Find everything within 2 Angstroms and say they are bonded. 94 else: 95 self._find_bonded_atoms(index, mol, radius=2) 96 97 # Loop over the bonded atoms. 98 matching_list = [] 99 for bonded_index in mol.bonded[index]: 100 if regex.search(mol.atom_name[bonded_index], attached_atom): 101 matching_list.append(bonded_index) 102 num_attached = len(matching_list) 103 104 # Problem. 105 if num_attached > 1: 106 # Get the atom names. 107 matching_names = [] 108 for i in matching_list: 109 matching_names.append(mol.atom_name[i]) 110 111 # Return nothing but a warning. 112 return None, None, None, None, None, 'More than one attached atom found: ' + repr(matching_names) 113 114 # No attached atoms. 115 if num_attached == 0: 116 return None, None, None, None, None, "No attached atom could be found" 117 118 # The bonded atom info. 119 index = matching_list[0] 120 bonded_num = mol.atom_num[index] 121 bonded_name = mol.atom_name[index] 122 element = mol.element[index] 123 pos = [mol.x[index], mol.y[index], mol.z[index]] 124 attached_name = mol.atom_name[index] 125 126 # Return the information. 127 return bonded_num, bonded_name, element, pos, attached_name, None
128 129
130 - def _find_bonded_atoms(self, index, mol, radius=1.2):
131 """Find all atoms within a sphere and say that they are attached to the central atom. 132 133 The found atoms will be added to the 'bonded' data structure. 134 135 136 @param index: The index of the central atom. 137 @type index: int 138 @param mol: The molecule container. 139 @type mol: MolContainer instance 140 """ 141 142 # Central atom info. 143 centre = array([mol.x[index], mol.y[index], mol.z[index]], float64) 144 145 # Atom loop. 146 dist_list = [] 147 connect_list = {} 148 element_list = {} 149 for i in range(len(mol.atom_num)): 150 # Skip proton to proton bonds! 151 if mol.element[index] == 'H' and mol.element[i] == 'H': 152 continue 153 154 # The atom's position. 155 pos = array([mol.x[i], mol.y[i], mol.z[i]], float64) 156 157 # The distance from the centre. 158 dist = linalg.norm(centre-pos) 159 160 # The atom is within the radius. 161 if dist < radius: 162 # Store the distance. 163 dist_list.append(dist) 164 165 # Store the atom index. 166 connect_list[dist] = i 167 168 # Store the element type. 169 element_list[dist] = mol.element[i] 170 171 # The maximum number of allowed covalent bonds. 172 max_conn = 1000 # Ridiculous default! 173 if mol.element[index] == 'H': 174 max_conn = 1 175 elif mol.element[index] == 'O': 176 max_conn = 2 177 elif mol.element[index] == 'N': 178 max_conn = 3 179 elif mol.element[index] == 'C': 180 max_conn = 4 181 182 # Sort. 183 dist_list.sort() 184 185 # Loop over the max number of connections (or the number of connected atoms, if less). 186 for i in range(min(max_conn, len(dist_list))): 187 mol.atom_connect(index, connect_list[dist_list[i]])
188 189
190 - def _get_chemical_name(self, hetID):
191 """Return the chemical name corresponding to the given residue ID. 192 193 The following names are currently returned:: 194 ________________________________________________ 195 | | | 196 | hetID | Chemical name | 197 |________|_____________________________________| 198 | | | 199 | TNS | Tensor | 200 | COM | Centre of mass | 201 | AXS | Tensor axes | 202 | SIM | Monte Carlo simulation tensor axes | 203 | PIV | Pivot point | 204 | CON | Cone object | 205 | AVE | Average vector | 206 |________|_____________________________________| 207 208 For any other residues, no description is returned. 209 210 @param hetID: The residue ID. 211 @type hetID: str 212 @return: The chemical name. 213 @rtype: str or None 214 """ 215 216 # Tensor. 217 if hetID == 'TNS': 218 return 'Tensor' 219 220 # Centre of mass. 221 if hetID == 'COM': 222 return 'Centre of mass' 223 224 # Tensor axes. 225 if hetID == 'AXS': 226 return 'Tensor axes' 227 228 # Monte Carlo simulation tensor axes. 229 if hetID == 'SIM': 230 return 'Monte Carlo simulation tensor axes' 231 232 # Pivot point. 233 if hetID == 'PIV': 234 return 'Pivot point' 235 236 # Cone object. 237 if hetID == 'CON': 238 return 'Cone' 239 240 # Average vector. 241 if hetID == 'AVE': 242 return 'Average vector'
243 244
245 - def _parse_models_gaussian(self, file_path, verbosity=1):
246 """Generator function for looping over the models in the Gaussian log file. 247 248 @param file_path: The full path of the Gaussian log file. 249 @type file_path: str 250 @return: The model number and all the records for that model. 251 @rtype: tuple of int and array of str 252 @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. 253 @type verbosity: int 254 """ 255 256 # Open the file. 257 file = open_read_file(file_path, verbosity=verbosity) 258 lines = file.readlines() 259 file.close() 260 261 # Check for empty files. 262 if lines == []: 263 raise RelaxError("The Gaussian log file is empty.") 264 265 # Init. 266 found = False 267 str_index = 0 268 total_atom = 0 269 model = 0 270 records = [] 271 272 # Loop over the data. 273 for i in range(len(lines)): 274 # Found a structure. 275 if search("Standard orientation", lines[i]): 276 found = True 277 str_index = 0 278 continue 279 280 # End of the model. 281 if found and str_index > 4 and search("---------", lines[i]): 282 # Yield the info 283 yield records 284 285 # Reset. 286 records = [] 287 found = False 288 289 # Not a structure line. 290 if not found: 291 continue 292 293 # Append the line as a record of the model. 294 records.append(lines[i]) 295 296 # Increment the structure line index. 297 str_index += 1
298 299
300 - def _parse_pdb_connectivity_annotation(self, lines):
301 """Loop over and parse the PDB connectivity annotation records. 302 303 These are the records identified in the U{PDB version 3.30 documentation<http://www.wwpdb.org/documentation/format33/sect6.html>}. 304 305 306 @param lines: The lines of the PDB file excluding the sections prior to the connectivity annotation section. 307 @type lines: list of str 308 @return: The remaining PDB lines with the connectivity annotation records stripped. 309 @rtype: list of str 310 """ 311 312 # The ordered list of record names in the connectivity annotation section. 313 records = [ 314 'SSBOND', 315 'LINK ', 316 'CISPEP' 317 ] 318 319 # Loop over the lines. 320 for i in range(len(lines)): 321 # No match, therefore assume to be out of the connectivity annotation section. 322 if lines[i][:6] not in records: 323 break 324 325 # Return the remaining lines. 326 return lines[i:]
327 328
329 - def _parse_pdb_coord(self, lines):
330 """Generator function for looping over the models in the PDB file. 331 332 These are the records identified in the PDB version 3.30 documentation at U{http://www.wwpdb.org/documentation/format33/sect9.html}. 333 334 335 @param lines: The lines of the coordinate section. 336 @type lines: list of str 337 @return: The model number and all the records for that model. 338 @rtype: tuple of int and array of str 339 """ 340 341 # Init. 342 model = None 343 records = [] 344 345 # Loop over the data. 346 for i in range(len(lines)): 347 # A new model record. 348 if lines[i][:5] == 'MODEL': 349 try: 350 model = int(lines[i].split()[1]) 351 except: 352 raise RelaxError("The MODEL record " + repr(lines[i]) + " is corrupt, cannot read the PDB file.") 353 354 # Skip all records prior to the first ATOM or HETATM record. 355 if not (lines[i][:4] == 'ATOM' or lines[i][:6] == 'HETATM') and not len(records): 356 continue 357 358 # End of the model. 359 if lines[i][:6] == 'ENDMDL': 360 # Yield the info. 361 yield model, records 362 363 # Reset the records. 364 records = [] 365 366 # Skip the rest of this loop. 367 continue 368 369 # Append the line as a record of the model. 370 records.append(lines[i]) 371 372 # If records is not empty then there are no models, so yield the lot. 373 if len(records): 374 yield model, records
375 376
377 - def _parse_pdb_hetrogen(self, lines):
378 """Loop over and parse the PDB hetrogen records. 379 380 These are the records identified in the PDB version 3.30 documentation at U{http://www.wwpdb.org/documentation/format33/sect4.html}. 381 382 383 @param lines: The lines of the PDB file excluding the sections prior to the hetrogen section. 384 @type lines: list of str 385 @return: The remaining PDB lines with the hetrogen records stripped. 386 @rtype: list of str 387 """ 388 389 # The ordered list of record names in the hetrogen section. 390 records = [ 391 'HET ', 392 'FORMUL', 393 'HETNAM', 394 'HETSYN' 395 ] 396 397 # Loop over the lines. 398 for i in range(len(lines)): 399 # No match, therefore assume to be out of the hetrogen section. 400 if lines[i][:6] not in records: 401 break 402 403 # Return the remaining lines. 404 return lines[i:]
405 406
407 - def _parse_pdb_misc(self, lines):
408 """Loop over and parse the PDB miscellaneous records. 409 410 These are the records identified in the PDB version 3.30 documentation at U{http://www.wwpdb.org/documentation/format33/sect7.html}. 411 412 413 @param lines: The lines of the PDB file excluding the sections prior to the miscellaneous section. 414 @type lines: list of str 415 @return: The remaining PDB lines with the miscellaneous records stripped. 416 @rtype: list of str 417 """ 418 419 # The ordered list of record names in the miscellaneous section. 420 records = [ 421 'SITE ' 422 ] 423 424 # Loop over the lines. 425 for i in range(len(lines)): 426 # No match, therefore assume to be out of the miscellaneous section. 427 if lines[i][:6] not in records: 428 break 429 430 # Return the remaining lines. 431 return lines[i:]
432 433
434 - def _parse_pdb_prim_struct(self, lines):
435 """Loop over and parse the PDB primary structure records. 436 437 These are the records identified in the PDB version 3.30 documentation at U{http://www.wwpdb.org/documentation/format33/sect3.html}. 438 439 440 @param lines: The lines of the PDB file excluding the title section. 441 @type lines: list of str 442 @return: The remaining PDB lines with the primary structure records stripped. 443 @rtype: list of str 444 """ 445 446 # The ordered list of record names in the primary structure section. 447 records = [ 448 'DBREF ', 449 'DBREF1', 450 'DBREF2', 451 'SEQADV', 452 'SEQRES', 453 'MODRES' 454 ] 455 456 # Loop over the lines. 457 for i in range(len(lines)): 458 # No match, therefore assume to be out of the primary structure section. 459 if lines[i][:6] not in records: 460 break 461 462 # Return the remaining lines. 463 return lines[i:]
464 465
466 - def _parse_pdb_ss(self, lines):
467 """Loop over and parse the PDB secondary structure records. 468 469 These are the records identified in the PDB version 3.30 documentation at U{http://www.wwpdb.org/documentation/format33/sect5.html}. 470 471 472 @param lines: The lines of the PDB file excluding the sections prior to the secondary structure section. 473 @type lines: list of str 474 @return: The remaining PDB lines with the secondary structure records stripped. 475 @rtype: list of str 476 """ 477 478 # The ordered list of record names in the secondary structure section (the depreciated TURN record is also included to handle old PDB files). 479 records = [ 480 'HELIX ', 481 'SHEET ', 482 'TURN ' 483 ] 484 485 # Loop over the lines. 486 for i in range(len(lines)): 487 # No match, therefore assume to be out of the secondary structure section. 488 if lines[i][:6] not in records: 489 break 490 491 # A helix. 492 if lines[i][:5] == 'HELIX': 493 # Parse the record. 494 record_type, ser_num, helix_id, init_res_name, init_chain_id, init_seq_num, init_icode, end_res_name, end_chain_id, end_seq_num, end_icode, helix_class, comment, length = pdb_read.helix(lines[i]) 495 496 # Store the data. 497 if not hasattr(self, 'helices'): 498 self.helices = [] 499 self.helices.append([helix_id, init_chain_id, init_res_name, init_seq_num, end_chain_id, end_res_name, end_seq_num, helix_class, length]) 500 501 # A sheet. 502 if lines[i][:5] == 'SHEET': 503 # Parse the record. 504 record_type, strand, sheet_id, num_strands, init_res_name, init_chain_id, init_seq_num, init_icode, end_res_name, end_chain_id, end_seq_num, end_icode, sense, cur_atom, cur_res_name, cur_chain_id, cur_res_seq, cur_icode, prev_atom, prev_res_name, prev_chain_id, prev_res_seq, prev_icode = pdb_read.sheet(lines[i]) 505 506 # Store the data. 507 if not hasattr(self, 'sheets'): 508 self.sheets = [] 509 self.sheets.append([strand, sheet_id, num_strands, init_res_name, init_chain_id, init_seq_num, init_icode, end_res_name, end_chain_id, end_seq_num, end_icode, sense, cur_atom, cur_res_name, cur_chain_id, cur_res_seq, cur_icode, prev_atom, prev_res_name, prev_chain_id, prev_res_seq, prev_icode]) 510 511 # Return the remaining lines. 512 return lines[i:]
513 514
515 - def _parse_pdb_title(self, lines):
516 """Loop over and parse the PDB title records. 517 518 These are the records identified in the PDB version 3.30 documentation at U{http://www.wwpdb.org/documentation/format33/sect2.html}. 519 520 521 @param lines: All lines of the PDB file. 522 @type lines: list of str 523 @return: The remaining PDB lines with the title records stripped. 524 @rtype: list of str 525 """ 526 527 # The ordered list of (sometimes truncated) record names in the title section. 528 records = [ 529 'HEADER', 530 'OBSLTE', 531 'TITLE ', 532 'SPLT ', 533 'CAVEAT', 534 'COMPND', 535 'SOURCE', 536 'KEYWDS', 537 'EXPDTA', 538 'NUMMDL', 539 'MDLTYP', 540 'AUTHOR', 541 'REVDAT', 542 'SPRSDE', 543 'JRNL ', 544 'REMARK' 545 ] 546 547 # Loop over the lines. 548 for i in range(len(lines)): 549 # No match, therefore assume to be out of the title section. 550 if lines[i][:6] not in records: 551 break 552 553 # Return the remaining lines. 554 return lines[i:]
555 556
557 - def _parse_pdb_transform(self, lines):
558 """Loop over and parse the PDB transform records. 559 560 These are the records identified in the PDB version 3.30 documentation at U{http://www.wwpdb.org/documentation/format33/sect8.html}. 561 562 563 @param lines: The lines of the PDB file excluding the sections prior to the transform section. 564 @type lines: list of str 565 @return: The remaining PDB lines with the transform records stripped. 566 @rtype: list of str 567 """ 568 569 # The ordered list of record names in the transform section. 570 records = [ 571 'CRYST', 572 'MTRIX', 573 'ORIGX', 574 'SCALE', 575 ] 576 577 # Loop over the lines. 578 for i in range(len(lines)): 579 # No match, therefore assume to be out of the transform section. 580 if lines[i][0: 5] not in records: 581 break 582 583 # Return the remaining lines. 584 return lines[i:]
585 586
587 - def _parse_models_xyz(self, file_path, verbosity=1):
588 """Generator function for looping over the models in the XYZ file. 589 590 @param file_path: The full path of the XYZ file. 591 @type file_path: str 592 @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. 593 @type verbosity: int 594 @return: The model number and all the records for that model. 595 @rtype: tuple of int and array of str 596 """ 597 598 # Open the file. 599 file = open_read_file(file_path, verbosity=verbosity) 600 lines = file.readlines() 601 file.close() 602 603 # Check for empty files. 604 if lines == []: 605 raise RelaxError("The XYZ file is empty.") 606 607 # Init. 608 total_atom = 0 609 model = 0 610 records = [] 611 612 # Loop over the data. 613 for i in range(len(lines)): 614 num=0 615 word = lines[i].split() 616 # Find the total atom number and the first model. 617 if (i==0) and (len(word)==1): 618 try: 619 total_atom = int(word[0]) 620 num = 1 621 except: 622 raise RelaxError("The MODEL record " + repr(lines[i]) + " is corrupt, cannot read the XYZ file.") 623 624 # End of the model. 625 if (len(records) == total_atom): 626 # Yield the info 627 yield records 628 629 # Reset the records. 630 records = [] 631 632 # Skip all records prior to atom coordinates record. 633 if (len(word) != 4): 634 continue 635 636 # Append the line as a record of the model. 637 records.append(lines[i]) 638 639 # If records is not empty then there are no models, so yield the lot. 640 if len(records): 641 yield records
642 643
644 - def _parse_mols(self, records):
645 """Generator function for looping over the molecules in the PDB records of a model. 646 647 @param records: The list of PDB records for the model, or if no models exist the entire PDB file. 648 @type records: list of str 649 @return: The molecule number and all the records for that molecule. 650 @rtype: tuple of int and list of str 651 """ 652 653 # Check for empty records. 654 if records == []: 655 raise RelaxError("There are no PDB records for this model.") 656 657 # Init. 658 mol_num = 1 659 mol_records = [] 660 end = False 661 662 # Loop over the data. 663 for i in range(len(records)): 664 # A PDB termination record. 665 if records[i][:3] == 'END': 666 break 667 668 # A master record, so we are done. 669 if records[i][:6] == 'MASTER': 670 break 671 672 # A model termination record. 673 if records[i][:6] == 'ENDMDL': 674 end = True 675 676 # A molecule termination record with no trailing HETATM. 677 elif i < len(records)-1 and records[i][:3] == 'TER' and not records[i+1][:6] == 'HETATM': 678 end = True 679 680 # A HETATM followed by an ATOM record. 681 elif i < len(records)-1 and records[i][:6] == 'HETATM' and records[i+1][:4] == 'ATOM': 682 end = True 683 684 # End. 685 if end: 686 # Yield the info. 687 yield mol_num, mol_records 688 689 # Reset the records. 690 mol_records = [] 691 692 # Increment the molecule number. 693 mol_num = mol_num + 1 694 695 # Reset the flag. 696 end = False 697 698 # Skip the rest of this loop. 699 continue 700 701 # Append the line as a record of the molecule. 702 mol_records.append(records[i]) 703 704 # If records is not empty then there is only a single molecule, so yield the lot. 705 if len(mol_records): 706 yield mol_num, mol_records
707 708
709 - def _pdb_chain_id_to_mol_index(self, chain_id=None):
710 """Convert the PDB chain ID into the molecule index in a regular way. 711 712 @keyword chain_id: The PDB chain ID string. 713 @type chain_id: str 714 @return: The corresponding molecule index. 715 @rtype: int 716 """ 717 718 # Initialise. 719 mol_index = 0 720 721 # Convert to the molecule index. 722 if chain_id: 723 mol_index = CHAIN_ID_LIST.index(chain_id) 724 725 # Return the index. 726 return mol_index
727 728
729 - def _residue_data(self, res_nums=None, res_names=None):
730 """Convert the residue info into a dictionary of unique residues with numbers as keys. 731 732 @keyword res_nums: The list of residue numbers. 733 @type res_nums: list of int 734 @keyword res_names: The list of residue names matching the numbers. 735 @type res_names: list of str 736 @return: The dictionary of residue names with residue numbers as keys. 737 @rtype: dict of str 738 """ 739 740 # Initialise. 741 data = {} 742 743 # Loop over the data. 744 for i in range(len(res_nums)): 745 # The residue data already exists. 746 if res_nums[i] in data: 747 continue 748 749 # Add the data. 750 data[res_nums[i]] = res_names[i] 751 752 # Return the dictionary. 753 return data
754 755
756 - def _validate_data_arrays(self, struct):
757 """Check the validity of the data arrays in the given structure object. 758 759 @param struct: The structural object. 760 @type struct: Structure_container instance 761 """ 762 763 # The number of atoms. 764 num = len(struct.atom_name) 765 766 # Check the other lengths. 767 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: 768 raise RelaxError("The structural data is invalid.")
769 770
771 - def _validate_records(self, lines):
772 """Make sure all PDB records are 80 char in length, padding with whitespace when needed. 773 774 All newline characters are stripped from the records as well. 775 776 777 @param lines: All lines of the PDB file. 778 @type lines: list of str 779 @return: The padded PDB lines. 780 @rtype: list of str 781 """ 782 783 # Loop over the lines. 784 for i in range(len(lines)): 785 # Strip the newline character. 786 lines[i] = lines[i].rstrip('\r\n') 787 788 # Pad if needed. 789 if len(lines[i]) != 80: 790 lines[i] = "%-80s" % lines[i] 791 792 # Return the fixed lines. 793 return lines
794 795
796 - def _mol_type(self, mol):
797 """Determine the type of molecule. 798 799 @param mol: The molecule data container. 800 @type mol: MolContainer instance 801 """ 802 803 # Amino acids. 804 aa = ['ALA', 'ARG', 'ASN', 'ASP', 'CYS', 'GLU', 'GLN', 'GLY', 'HIS', 'ILE', 'LEU', 'LYS', 'MET', 'PHE', 'PRO', 'SER', 'THR', 'TRP', 'TYR', 'VAL'] 805 806 # Set the molecule type to default to 'other'. 807 mol.type = 'other' 808 809 # Loop over the residues. 810 for res in mol.res_name: 811 # Protein. 812 if res in aa: 813 # Set the molecule type and return. 814 mol.type = 'protein' 815 return
816 817
818 - def _protein_connect(self, mol):
819 """Set up the connectivities for the protein. 820 821 @param mol: The molecule data container. 822 @type mol: MolContainer instance 823 """ 824 825 # Initialise some residue data. 826 curr_res_num = None 827 res_atoms = [] 828 829 # Loop over all atoms. 830 for i in range(len(mol.atom_num)): 831 # New residue. 832 if mol.res_num[i] != curr_res_num: 833 # Intra-residue connectivites. 834 if len(res_atoms): 835 self._protein_intra_connect(mol, res_atoms) 836 837 # Update the residue number. 838 curr_res_num = mol.res_num[i] 839 840 # Reset the residue atom index list. 841 res_atoms = [] 842 843 # Add the atom index to the list. 844 res_atoms.append(i) 845 846 # Last atom. 847 if i == len(mol.atom_num) - 1 and len(res_atoms): 848 self._protein_intra_connect(mol, res_atoms)
849 850
851 - def _protein_intra_connect(self, mol, res_atoms):
852 """Set up the connectivities for the protein. 853 854 @param mol: The molecule data container. 855 @type mol: MolContainer instance 856 @param res_atoms: The list of atom indices corresponding to the residue. 857 @type res_atoms: list of int 858 """ 859 860 # Back bond connectivity. 861 indices = { 862 'N': None, 863 'C': None, 864 'O': None, 865 'CA': None, 866 'HN': None, 867 'H': None, # Same as HN. 868 'HA': None 869 } 870 871 # Loop over all atoms to find the indices. 872 for index in res_atoms: 873 if mol.atom_name[index] in indices: 874 indices[mol.atom_name[index]] = index 875 876 # Connect the atom pairs. 877 pairs = [ 878 ['N', 'HN'], 879 ['N', 'H'], 880 ['N', 'CA'], 881 ['CA', 'HA'], 882 ['CA', 'C'], 883 ['C', 'O'] 884 ] 885 886 # Loop over the atoms pairs and connect them. 887 for pair in pairs: 888 if indices[pair[0]] != None and indices[pair[1]] != None: 889 mol.atom_connect(indices[pair[0]], indices[pair[1]])
890 891
892 - def _translate(self, data, format='str'):
893 """Convert the data into a format for writing to file. 894 895 @param data: The data to convert to the required format. 896 @type data: anything 897 @keyword format: The format to convert to. This can be 'str', 'float', or 'int'. 898 @type format: str 899 @return: The converted version of the data. 900 @rtype: str 901 """ 902 903 # Conversion to string. 904 if format == 'str': 905 # None values. 906 if data == None: 907 data = '' 908 909 # Force convert to string. 910 if not isinstance(data, str): 911 data = repr(data) 912 913 # Conversion to float. 914 if format == 'float': 915 # None values. 916 if data == None: 917 data = 0.0 918 919 # Force convert to float. 920 if not isinstance(data, float): 921 data = float(data) 922 923 # Return the converted data. 924 return data
925 926
927 - def _trim_helix(self, helix=None, trim_res_list=[], res_data=None):
928 """Trim the given helix based on the list of deleted residue numbers. 929 930 @keyword helix: The single helix metadata structure. 931 @type helix: list 932 @keyword trim_res_list: The list of residue numbers which no longer exist. 933 @type trim_res_list: list of int 934 @keyword res_data: The dictionary of residue names with residue numbers as keys. 935 @type res_data: dict of str 936 @return: The trimmed helix metadata structure, or None if the whole helix is to be deleted. 937 @rtype: list or None 938 """ 939 940 # Unpack the helix residue numbers. 941 start_res = helix[3] 942 end_res = helix[6] 943 944 # The reverse residue list. 945 trim_res_list_rev = deepcopy(trim_res_list) 946 trim_res_list_rev.reverse() 947 948 # The helix residues. 949 helix_res = list(range(start_res, end_res+1)) 950 951 # Trim forwards. 952 for res_num in trim_res_list: 953 if res_num == start_res: 954 # Remove the residue. 955 helix_res.pop(0) 956 957 # No helix left. 958 if len(helix_res) == 0: 959 break 960 961 # Realias the starting residue. 962 start_res = helix_res[0] 963 964 # No helix left. 965 if len(helix_res) == 0: 966 return None 967 968 # Trim backwards. 969 for res_num in trim_res_list_rev: 970 if res_num == end_res: 971 helix_res.pop(-1) 972 end_res = helix_res[-1] 973 974 # Replace the starting and ending residues. 975 if start_res != helix[3]: 976 helix[3] = start_res 977 helix[2] = res_data[start_res] 978 if end_res != helix[6]: 979 helix[6] = end_res 980 helix[5] = res_data[end_res] 981 982 # The helix length. 983 helix[-1] = len(helix_res) 984 985 # Return the modified helix. 986 return helix
987 988
989 - def _trim_sheet(self, sheet=None, trim_res_list=[], res_data=None):
990 """Trim the given sheet based on the list of deleted residue numbers. 991 992 @keyword sheet: The single sheet metadata structure. 993 @type sheet: list 994 @keyword trim_res_list: The list of residue numbers which no longer exist. 995 @type trim_res_list: list of int 996 @keyword res_data: The dictionary of residue names with residue numbers as keys. 997 @type res_data: dict of str 998 @return: The trimmed sheet metadata structure, or None if the whole sheet is to be deleted. 999 @rtype: list or None 1000 """ 1001 1002 # Unpack the sheet residue numbers. 1003 start_res = sheet[5] 1004 end_res = sheet[9] 1005 1006 # The reverse residue list. 1007 trim_res_list_rev = deepcopy(trim_res_list) 1008 trim_res_list_rev.reverse() 1009 1010 # The sheet residues. 1011 sheet_res = list(range(start_res, end_res+1)) 1012 1013 # Trim forwards. 1014 for res_num in trim_res_list: 1015 if res_num == start_res: 1016 # Remove the residue. 1017 sheet_res.pop(0) 1018 1019 # No sheet left. 1020 if len(sheet_res) == 0: 1021 break 1022 1023 # Realias the starting residue. 1024 start_res = sheet_res[0] 1025 1026 # No sheet left. 1027 if len(sheet_res) == 0: 1028 return None 1029 1030 # Trim backwards. 1031 for res_num in trim_res_list_rev: 1032 if res_num == end_res: 1033 sheet_res.pop(-1) 1034 end_res = sheet_res[-1] 1035 1036 # Replace the starting and ending residues. 1037 if start_res != sheet[5]: 1038 sheet[5] = start_res 1039 sheet[3] = res_data[start_res] 1040 if end_res != sheet[9]: 1041 sheet[9] = end_res 1042 sheet[7] = res_data[end_res] 1043 1044 # Return the modified sheet. 1045 return sheet
1046 1047
1048 - 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):
1049 """Add a new atom to the structural data object. 1050 1051 @keyword mol_name: The name of the molecule. 1052 @type mol_name: str 1053 @keyword atom_name: The atom name, e.g. 'H1'. 1054 @type atom_name: str or None 1055 @keyword res_name: The residue name. 1056 @type res_name: str or None 1057 @keyword res_num: The residue number. 1058 @type res_num: int or None 1059 @keyword pos: The position vector of coordinates. If a rank-2 array is supplied, the length of the first dimension must match the number of models. 1060 @type pos: rank-1 or rank-2 array or list of float 1061 @keyword element: The element symbol. 1062 @type element: str or None 1063 @keyword atom_num: The atom number. 1064 @type atom_num: int or None 1065 @keyword chain_id: The chain identifier. 1066 @type chain_id: str or None 1067 @keyword segment_id: The segment identifier. 1068 @type segment_id: str or None 1069 @keyword pdb_record: The optional PDB record name, e.g. 'ATOM' or 'HETATM'. 1070 @type pdb_record: str or None 1071 """ 1072 1073 # Add a model if not present. 1074 if len(self.structural_data) == 0: 1075 self.add_model() 1076 1077 # Check the position. 1078 if is_float(pos[0]): 1079 if len(pos) != 3: 1080 raise RelaxError("The single atomic position %s must be a 3D list." % pos) 1081 else: 1082 if len(pos) != len(self.structural_data): 1083 raise RelaxError("The %s atomic positions does not match the %s models present." % (len(pos), len(self.structural_data))) 1084 1085 # Loop over each model. 1086 for i in range(len(self.structural_data)): 1087 # Alias the model. 1088 model = self.structural_data[i] 1089 1090 # Specific molecule. 1091 mol = self.get_molecule(mol_name, model=model.num) 1092 1093 # Add the molecule, if it does not exist. 1094 if mol == None: 1095 self.add_molecule(name=mol_name) 1096 mol = self.get_molecule(mol_name, model=model.num) 1097 1098 # Split up the position if needed. 1099 if is_float(pos[0]): 1100 model_pos = pos 1101 else: 1102 model_pos = pos[i] 1103 1104 # Add the atom. 1105 mol.atom_add(atom_name=atom_name, res_name=res_name, res_num=res_num, pos=model_pos, element=element, atom_num=atom_num, chain_id=chain_id, segment_id=segment_id, pdb_record=pdb_record)
1106 1107
1108 - def add_model(self, model=None, coords_from=None):
1109 """Add a new model to the store. 1110 1111 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. 1112 1113 @keyword model: The number of the model to create. 1114 @type model: int or None 1115 @keyword coords_from: The model number to take the coordinates from. 1116 @type coords_from: int or None 1117 @return: The model container. 1118 @rtype: ModelContainer instance 1119 """ 1120 1121 # Check if the model currently exists. 1122 if model != None: 1123 for i in range(len(self.structural_data)): 1124 if model == self.structural_data[i].num: 1125 raise RelaxError("The model '%s' already exists." % model) 1126 1127 # Add a new model. 1128 self.structural_data.add_item(model_num=model) 1129 1130 # The model to duplicate. 1131 if coords_from == None: 1132 coords_from = self.structural_data[0].num 1133 1134 # Construct the structural data for the model from the other models. 1135 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): 1136 # Add the atom. 1137 self.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) 1138 1139 # Return the model. 1140 return self.structural_data[-1]
1141 1142
1143 - def add_molecule(self, model_num=None, name=None):
1144 """Add a new molecule to the store. 1145 1146 @keyword model_num: The optional model to add the molecule to. If not supplied, the molecule will be added to all models. 1147 @type model_num: None or int 1148 @keyword name: The molecule identifier string. 1149 @type name: str 1150 """ 1151 1152 # Add a model if necessary. 1153 if len(self.structural_data) == 0: 1154 self.add_model() 1155 1156 # Add the molecule to each model. 1157 for model in self.model_loop(model=model_num): 1158 model.mol.add_item(mol_name=name, mol_cont=MolContainer())
1159 1160
1161 - def are_bonded(self, atom_id1=None, atom_id2=None):
1162 """Determine if two atoms are directly bonded to each other. 1163 1164 @keyword atom_id1: The molecule, residue, and atom identifier string of the first atom. 1165 @type atom_id1: str 1166 @keyword atom_id2: The molecule, residue, and atom identifier string of the second atom. 1167 @type atom_id2: str 1168 @return: True if the atoms are directly bonded. 1169 @rtype: bool 1170 """ 1171 1172 # Generate the selection objects. 1173 sel_obj1 = Selection(atom_id1) 1174 sel_obj2 = Selection(atom_id2) 1175 1176 # Build the connectivities if needed. 1177 for mol in self.structural_data[0].mol: 1178 for i in range(len(mol.atom_num)): 1179 if not len(mol.bonded[i]): 1180 self._find_bonded_atoms(i, mol, radius=2) 1181 1182 # Loop over the molecules. 1183 for mol in self.structural_data[0].mol: 1184 # Skip non-matching molecules. 1185 if not sel_obj1.contains_mol(mol.mol_name): 1186 continue 1187 if not sel_obj2.contains_mol(mol.mol_name): 1188 continue 1189 1190 # Find the first atom. 1191 index1 = None 1192 index2 = None 1193 for i in range(len(mol.atom_num)): 1194 # Matching first atom. 1195 if index1 == None and sel_obj1.contains_spin(mol.atom_num[i], mol.atom_name[i], mol.res_num[i], mol.res_name[i], mol.mol_name): 1196 index1 = i 1197 1198 # Matching second atom. 1199 if index2 == None and sel_obj2.contains_spin(mol.atom_num[i], mol.atom_name[i], mol.res_num[i], mol.res_name[i], mol.mol_name): 1200 index2 = i 1201 1202 # Nothing left to do. 1203 if index1 != None and index2 != None: 1204 break 1205 1206 # Connectivities exist. 1207 if index1 < len(mol.bonded): 1208 if index2 in mol.bonded[index1]: 1209 return True 1210 else: 1211 return False
1212 1213
1214 - def are_bonded_index(self, mol_index1=None, atom_index1=None, mol_index2=None, atom_index2=None):
1215 """Determine if two atoms, given as indices, are directly bonded to each other. 1216 1217 @keyword mol_index1: The molecule index of the first atom. 1218 @type mol_index1: int 1219 @keyword atom_index1: The index of the first atom. 1220 @type atom_index1: int 1221 @keyword mol_index2: The molecule index of the second atom. 1222 @type mol_index2: int 1223 @keyword atom_index2: The index of the second atom. 1224 @type atom_index2: int 1225 @return: True if the atoms are directly bonded. 1226 @rtype: bool 1227 """ 1228 1229 # Alias the molecule. 1230 mol1 = self.structural_data[0].mol[mol_index1] 1231 mol2 = self.structural_data[0].mol[mol_index2] 1232 1233 # Build the connectivities if needed. 1234 if not len(mol1.bonded[atom_index1]): 1235 self._find_bonded_atoms(atom_index1, mol1, radius=2) 1236 if not len(mol2.bonded[atom_index2]): 1237 self._find_bonded_atoms(atom_index2, mol2, radius=2) 1238 1239 # Is the second atom in the bonded list of the first? 1240 if atom_index2 in mol1.bonded[atom_index1]: 1241 return True 1242 1243 # No! 1244 return False
1245 1246
1247 - def atom_loop(self, atom_id=None, str_id=None, model_num=None, 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, mol_index_flag=False, index_flag=False, ave=False):
1248 """Generator function for looping over all atoms in the internal relax structural object. 1249 1250 This method should be designed as a U{generator<http://www.python.org/dev/peps/pep-0255/>}. It should loop over all atoms of the system yielding the following atomic information, if the corresponding flag is True, in tuple form: 1251 1252 1. Model number. 1253 2. Molecule name. 1254 3. Residue number. 1255 4. Residue name. 1256 5. Atom number. 1257 6. Atom name. 1258 7. The element name (its atomic symbol and optionally the isotope, e.g. 'N', 'Mg', '17O', '13C', etc). 1259 8. The position of the atom in Euclidean space. 1260 1261 1262 @keyword atom_id: The molecule, residue, and atom identifier string. Only atoms matching this selection will be yielded. 1263 @type atom_id: str 1264 @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. 1265 @type str_id: str, int, or None 1266 @keyword model_num: Only loop over a specific model. 1267 @type model_num: int or None 1268 @keyword mol_name_flag: A flag which if True will cause the molecule name to be yielded. 1269 @type mol_name_flag: bool 1270 @keyword res_num_flag: A flag which if True will cause the residue number to be yielded. 1271 @type res_num_flag: bool 1272 @keyword res_name_flag: A flag which if True will cause the residue name to be yielded. 1273 @type res_name_flag: bool 1274 @keyword atom_num_flag: A flag which if True will cause the atom number to be yielded. 1275 @type atom_num_flag: bool 1276 @keyword atom_name_flag: A flag which if True will cause the atom name to be yielded. 1277 @type atom_name_flag: bool 1278 @keyword element_flag: A flag which if True will cause the element name to be yielded. 1279 @type element_flag: bool 1280 @keyword pos_flag: A flag which if True will cause the atomic position to be yielded. 1281 @type pos_flag: bool 1282 @keyword mol_index_flag: A flag which if True will cause the molecule index to be yielded. 1283 @type mol_index_flag: bool 1284 @keyword index_flag: A flag which if True will cause the atomic index to be yielded. 1285 @type index_flag: bool 1286 @keyword ave: A flag which if True will result in this method returning the average atom properties across all loaded structures. 1287 @type ave: bool 1288 @return: A tuple of atomic information, as described in the docstring. 1289 @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). 1290 """ 1291 1292 # Check that the structure is loaded. 1293 if not len(self.structural_data): 1294 raise RelaxNoPdbError 1295 1296 # Generate the selection object. 1297 sel_obj = None 1298 if atom_id: 1299 sel_obj = Selection(atom_id) 1300 1301 # Obtain all data from the first model (except the position data). 1302 model = self.structural_data[0] 1303 1304 # Loop over the molecules. 1305 for mol_index in range(len(model.mol)): 1306 mol = model.mol[mol_index] 1307 1308 # Skip non-matching molecules. 1309 if sel_obj and not sel_obj.contains_mol(mol.mol_name): 1310 continue 1311 1312 # Loop over all atoms. 1313 for i in range(len(mol.atom_name)): 1314 # Skip non-matching atoms. 1315 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): 1316 continue 1317 1318 # Initialise. 1319 res_num = mol.res_num[i] 1320 res_name = mol.res_name[i] 1321 atom_num = mol.atom_num[i] 1322 atom_name = mol.atom_name[i] 1323 element = mol.element[i] 1324 1325 # The atom position. 1326 if pos_flag: 1327 # Average the position. 1328 if ave: 1329 # Initialise. 1330 pos = zeros(3, float64) 1331 1332 # Loop over the models. 1333 for model in self.model_loop(model=model_num): 1334 # Alias. 1335 mol2 = model.mol[mol_index] 1336 1337 # Some sanity checks. 1338 if mol2.atom_num[i] != atom_num: 1339 raise RelaxError("The loaded structures do not contain the same atoms. The average structural properties can not be calculated.") 1340 1341 # Sum the atom positions. 1342 pos = pos + array([mol2.x[i], mol2.y[i], mol2.z[i]], float64) 1343 1344 # Average the position array (divide by the number of models). 1345 pos = pos / len(self.structural_data) 1346 1347 # All positions. 1348 else: 1349 # Initialise. 1350 pos = [] 1351 1352 # Loop over the models. 1353 for model in self.model_loop(model=model_num): 1354 # Alias. 1355 mol2 = model.mol[mol_index] 1356 1357 # Append the position. 1358 pos.append([mol2.x[i], mol2.y[i], mol2.z[i]]) 1359 1360 # Convert. 1361 pos = array(pos, float64) 1362 1363 # The molecule name. 1364 mol_name = mol.mol_name 1365 1366 # Build the tuple to be yielded. 1367 atomic_tuple = () 1368 if mol_name_flag: 1369 atomic_tuple = atomic_tuple + (mol_name,) 1370 if res_num_flag: 1371 atomic_tuple = atomic_tuple + (res_num,) 1372 if res_name_flag: 1373 atomic_tuple = atomic_tuple + (res_name,) 1374 if atom_num_flag: 1375 atomic_tuple = atomic_tuple + (atom_num,) 1376 if atom_name_flag: 1377 atomic_tuple = atomic_tuple + (atom_name,) 1378 if element_flag: 1379 atomic_tuple = atomic_tuple + (element,) 1380 if pos_flag: 1381 atomic_tuple = atomic_tuple + (pos,) 1382 if mol_index_flag: 1383 atomic_tuple += (mol_index,) 1384 if index_flag: 1385 atomic_tuple += (i,) 1386 1387 # Yield the information. 1388 if len(atomic_tuple) == 1: 1389 atomic_tuple = atomic_tuple[0] 1390 yield atomic_tuple
1391 1392
1393 - 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):
1394 """Find the bond vector between the atoms of 'attached_atom' and 'atom_id'. 1395 1396 @keyword attached_atom: The name of the bonded atom. 1397 @type attached_atom: str 1398 @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. 1399 @type model_num: None or int 1400 @keyword mol_name: The name of the molecule that attached_atom belongs to. 1401 @type mol_name: str 1402 @keyword res_num: The number of the residue that attached_atom belongs to. 1403 @type res_num: str 1404 @keyword res_name: The name of the residue that attached_atom belongs to. 1405 @type res_name: str 1406 @keyword spin_num: The number of the spin that attached_atom is attached to. 1407 @type spin_num: str 1408 @keyword spin_name: The name of the spin that attached_atom is attached to. 1409 @type spin_name: str 1410 @keyword return_name: A flag which if True will cause the name of the attached atom to be returned together with the bond vectors. 1411 @type return_name: bool 1412 @keyword return_warnings: A flag which if True will cause warning messages to be returned. 1413 @type return_warnings: bool 1414 @return: The list of bond vectors for each model. 1415 @rtype: list of numpy arrays (or a tuple if return_name or return_warnings are set) 1416 """ 1417 1418 # Initialise some objects. 1419 vectors = [] 1420 attached_name = None 1421 warnings = None 1422 1423 # Use the first model for the atom matching. 1424 model = self.structural_data[0] 1425 1426 # Loop over the molecules. 1427 for mol_index in range(len(model.mol)): 1428 # Alias. 1429 mol = model.mol[mol_index] 1430 1431 # Skip non-matching molecules. 1432 if mol_name and mol_name != mol.mol_name: 1433 continue 1434 1435 # Find the atomic index of the base atom. 1436 index = None 1437 for i in range(len(mol.atom_name)): 1438 # Residues don't match. 1439 if (res_num != None and mol.res_num[i] != res_num) or (res_name != None and mol.res_name[i] != res_name): 1440 continue 1441 1442 # Atoms don't match. 1443 if (spin_num != None and mol.atom_num[i] != spin_num) or (spin_name != None and mol.atom_name[i] != spin_name): 1444 continue 1445 1446 # Update the index and stop searching. 1447 index = i 1448 break 1449 1450 # Found the atom. 1451 if index != None: 1452 # Loop over the models. 1453 for model in self.model_loop(model=model_num): 1454 # Alias. 1455 mol = model.mol[mol_index] 1456 1457 # Get the atom bonded to this model/molecule/residue/atom. 1458 bonded_num, bonded_name, element, pos, attached_name, warnings = self._bonded_atom(attached_atom, index, mol) 1459 1460 # No bonded atom. 1461 if (bonded_num, bonded_name, element) == (None, None, None): 1462 continue 1463 1464 # The bond vector. 1465 vector = array(pos, float64) - array([mol.x[index], mol.y[index], mol.z[index]], float64) 1466 1467 # Append the vector to the vectors array. 1468 vectors.append(vector) 1469 1470 # Not found. 1471 else: 1472 warnings = "Cannot find the atom in the structure" 1473 1474 # Build the tuple to be yielded. 1475 data = (vectors,) 1476 if return_name: 1477 data = data + (attached_name,) 1478 if return_warnings: 1479 data = data + (warnings,) 1480 1481 # Return the data. 1482 return data
1483 1484
1485 - def connect_atom(self, mol_name=None, index1=None, index2=None):
1486 """Connect two atoms in the structural data object. 1487 1488 @keyword mol_name: The name of the molecule. 1489 @type mol_name: str 1490 @keyword index1: The global index of the first atom. 1491 @type index1: str 1492 @keyword index2: The global index of the first atom. 1493 @type index2: str 1494 """ 1495 1496 # Add the molecule, if it does not exist. 1497 if self.get_molecule(mol_name) == None: 1498 self.add_molecule(name=mol_name) 1499 1500 # Loop over each model. 1501 for model in self.structural_data: 1502 # Specific molecule. 1503 mol = self.get_molecule(mol_name) 1504 1505 # Add the atom. 1506 mol.atom_connect(index1=index1, index2=index2)
1507 1508
1509 - def delete(self, atom_id=None, verbosity=1):
1510 """Deletion of structural information. 1511 1512 @keyword atom_id: The molecule, residue, and atom identifier string. This matches the spin ID string format. If not given, then all structural data will be deleted. 1513 @type atom_id: str or None 1514 @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. 1515 @type verbosity: int 1516 """ 1517 1518 # All data. 1519 if atom_id == None: 1520 # Printout. 1521 if verbosity: 1522 print("Deleting the following structural data:\n") 1523 print(self.structural_data) 1524 1525 # Delete the structural data. 1526 del self.structural_data 1527 1528 # Initialise the empty model list. 1529 self.structural_data = ModelList() 1530 1531 # Atom subset deletion. 1532 else: 1533 # Generate the selection object. 1534 sel_obj = None 1535 if atom_id: 1536 sel_obj = Selection(atom_id) 1537 1538 # Loop over the atoms and find the indices of the atoms to delete. 1539 indices = [] 1540 for i in self.atom_loop(atom_id=atom_id, index_flag=True): 1541 indices.append(i) 1542 indices.reverse() 1543 1544 # Loop over the models. 1545 del_res_nums = [] 1546 for model in self.model_loop(): 1547 # Loop over the molecules. 1548 for mol_index in range(len(model.mol)): 1549 mol = model.mol[mol_index] 1550 1551 # Skip non-matching molecules. 1552 if sel_obj and not sel_obj.contains_mol(mol.mol_name): 1553 continue 1554 1555 # Generate a residue data dictionary for the metadata trimming (prior to atom deletion). 1556 res_data = self._residue_data(res_nums=mol.res_num, res_names=mol.res_name) 1557 1558 # Loop over the reverse indices and pop out the data. 1559 for i in indices: 1560 mol.atom_num.pop(i) 1561 mol.atom_name.pop(i) 1562 mol.bonded.pop(i) 1563 mol.chain_id.pop(i) 1564 mol.element.pop(i) 1565 mol.pdb_record.pop(i) 1566 mol.res_name.pop(i) 1567 res_num = mol.res_num.pop(i) 1568 mol.seg_id.pop(i) 1569 mol.x.pop(i) 1570 mol.y.pop(i) 1571 mol.z.pop(i) 1572 1573 # The residue no longer exists. 1574 if res_num not in mol.res_num and res_num not in del_res_nums: 1575 del_res_nums.append(res_num) 1576 1577 # Nothing more to do. 1578 if not len(del_res_nums): 1579 return 1580 1581 # Fix the deleted residue number order. 1582 del_res_nums.reverse() 1583 1584 # Handle the helix metadata. 1585 if hasattr(self, 'helices'): 1586 del_helix_indices = [] 1587 for i in range(len(self.helices)): 1588 # Trim the helix. 1589 helix = self._trim_helix(helix=self.helices[i], trim_res_list=del_res_nums, res_data=res_data) 1590 1591 # Trimmed helix. 1592 if helix != None: 1593 self.helices[i] = helix 1594 1595 # No helix left. 1596 else: 1597 del_helix_indices.append(i) 1598 1599 # Loop over the reverse helix indices and pop out the data. 1600 del_helix_indices.reverse() 1601 for i in del_helix_indices: 1602 self.helices.pop(i) 1603 1604 # Handle the sheet metadata. 1605 if hasattr(self, 'sheets'): 1606 del_sheet_indices = [] 1607 for i in range(len(self.sheets)): 1608 # Trim the sheet. 1609 sheet = self._trim_sheet(sheet=self.sheets[i], trim_res_list=del_res_nums, res_data=res_data) 1610 1611 # Trimmed sheet. 1612 if sheet != None: 1613 self.sheets[i] = sheet 1614 1615 # No sheet left. 1616 else: 1617 del_sheet_indices.append(i) 1618 1619 # Loop over the reverse sheet indices and pop out the data. 1620 del_sheet_indices.reverse() 1621 for i in del_sheet_indices: 1622 self.sheets.pop(i)
1623 1624
1625 - def empty(self):
1626 """Report if the structural data structure is empty or not. 1627 1628 @return: True if empty, False otherwise. 1629 @rtype: bool 1630 """ 1631 1632 # Check the ModelList structure. 1633 if len(self.structural_data) == 0: 1634 return True 1635 else: 1636 return False
1637 1638
1639 - def from_xml(self, str_node, dir=None, file_version=1):
1640 """Recreate the structural object from the XML structural object node. 1641 1642 @param str_node: The structural object XML node. 1643 @type str_node: xml.dom.minicompat.Element instance 1644 @keyword dir: The name of the directory containing the results file. 1645 @type dir: str 1646 @keyword file_version: The relax XML version of the XML file. 1647 @type file_version: int 1648 """ 1649 1650 # Recreate all base objects (i.e. metadata). 1651 xml_to_object(str_node, self, file_version=file_version, blacklist=['model', 'displacements']) 1652 1653 # Recreate the model / molecule data structure. 1654 model_nodes = str_node.getElementsByTagName('model') 1655 self.structural_data.from_xml(model_nodes, file_version=file_version) 1656 1657 # The displacement structure. 1658 disp_nodes = str_node.getElementsByTagName('displacements') 1659 if len(disp_nodes): 1660 # Initialise the object. 1661 self.displacements = Displacements() 1662 1663 # Recreate the molecule data structures for the current model. 1664 self.displacements.from_xml(disp_nodes[0], file_version=file_version)
1665 1666
1667 - def get_model(self, model):
1668 """Return or create the model. 1669 1670 @param model: The model number. 1671 @type model: int or None 1672 @return: The ModelContainer corresponding to the model number or that newly created. 1673 @rtype: ModelContainer instance 1674 """ 1675 1676 # Check if the target is a single model. 1677 if model == None and self.num_models() > 1: 1678 raise RelaxError("The target model cannot be determined as there are %s models already present." % self.num_modes()) 1679 1680 # No model specified. 1681 if model == None: 1682 # Create the first model, if necessary. 1683 if self.num_models(): 1684 self.structural_data.add_item() 1685 1686 # Alias the first model. 1687 model_cont = self.structural_data[0] 1688 1689 # The model has been specified. 1690 else: 1691 # Get the preexisting model. 1692 found = False 1693 for model_cont in self.structural_data: 1694 if model_cont.num == model: 1695 found = True 1696 break 1697 1698 # Add the model if it doesn't exist. 1699 if not found: 1700 self.structural_data.add_item(model) 1701 model_cont = self.structural_data[-1]
1702 1703
1704 - def get_molecule(self, molecule, model=None):
1705 """Return the molecule. 1706 1707 Only one model can be specified. 1708 1709 1710 @param molecule: The molecule name. 1711 @type molecule: int or None 1712 @keyword model: The model number. 1713 @type model: int or None 1714 @raises RelaxError: If the model is not specified and there is more than one model loaded. 1715 @return: The MolContainer corresponding to the molecule name and model number. 1716 @rtype: MolContainer instance or None 1717 """ 1718 1719 # Check if the target is a single molecule. 1720 if model == None and self.num_models() > 1: 1721 raise RelaxError("The target molecule cannot be determined as there are %s models already present." % self.num_models()) 1722 1723 # Check the model argument. 1724 if not isinstance(model, int) and not model == None: 1725 raise RelaxNoneIntError 1726 1727 # No models. 1728 if not len(self.structural_data): 1729 return 1730 1731 # Loop over the models. 1732 for model_cont in self.model_loop(model): 1733 # Loop over the molecules. 1734 for mol in model_cont.mol: 1735 # Return the matching molecule. 1736 if mol.mol_name == molecule: 1737 return mol
1738 1739
1740 - def has_molecule(self, model_num=None, name=None):
1741 """Check if the molecule name exists. 1742 1743 @keyword model_num: The optional model to check. If not supplied, the molecule will be searched for across all models. 1744 @type model_num: None or int 1745 @param name: The molecule name. 1746 @type name: str 1747 @return: True if the molecule exists, False otherwise. 1748 @rtype: bool 1749 """ 1750 1751 # No models. 1752 if not len(self.structural_data): 1753 return False 1754 1755 # Loop over the models. 1756 for model_cont in self.model_loop(model=model_num): 1757 # Loop over the molecules. 1758 for mol in model_cont.mol: 1759 # Matching molecule. 1760 if mol.mol_name == name: 1761 return True 1762 1763 # No match. 1764 return False
1765 1766
1767 - def load_gaussian(self, file_path, set_mol_name=None, set_model_num=None, verbosity=False):
1768 """Method for loading structures from a Gaussian log file. 1769 1770 @param file_path: The full path of the Gaussian log file. 1771 @type file_path: str 1772 @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. 1773 @type set_mol_name: None, str, or list of str 1774 @keyword set_model_num: Set the model number of the loaded molecule. If set to None, then the Gaussian log model numbers will be preserved, if they exist. 1775 @type set_model_num: None, int, or list of int 1776 @keyword verbosity: A flag which if True will cause messages to be printed. 1777 @type verbosity: bool 1778 @return: The status of the loading of the Gaussian log file. 1779 @rtype: bool 1780 """ 1781 1782 # Initial printout. 1783 if verbosity: 1784 print("\nInternal relax Gaussian log parser.") 1785 1786 # Test if the file exists. 1787 if not access(file_path, F_OK): 1788 # Exit indicating failure. 1789 return False 1790 1791 # Separate the file name and path. 1792 path, file_name = os.path.split(file_path) 1793 1794 # The absolute path. 1795 path_abs = abspath(curdir) + sep + path 1796 1797 # Set up the molecule name data structure. 1798 if set_mol_name: 1799 if not isinstance(set_mol_name, list): 1800 set_mol_name = [set_mol_name] 1801 else: 1802 set_mol_name = [file_root(file_name) + '_mol1'] 1803 1804 # Set up the model number data structure. 1805 if set_model_num: 1806 if not isinstance(set_model_num, list): 1807 set_model_num = [set_model_num] 1808 else: 1809 set_model_num = [1] 1810 1811 # Loop over all models in the Gaussian log file, doing nothing so the last model records are stored. 1812 for model_records in self._parse_models_gaussian(file_path, verbosity=verbosity): 1813 pass 1814 1815 # Generate the molecule container. 1816 mol = MolContainer() 1817 1818 # Fill the molecular data object. 1819 mol.fill_object_from_gaussian(model_records) 1820 1821 # Create the structural data data structures. 1822 self.pack_structs([[mol]], orig_model_num=[1], set_model_num=set_model_num, orig_mol_num=[0], set_mol_name=set_mol_name, file_name=file_name, file_path=path, file_path_abs=path_abs, verbosity=verbosity) 1823 1824 # Loading worked. 1825 return True
1826 1827
1828 - def load_pdb(self, file_path, read_mol=None, set_mol_name=None, read_model=None, set_model_num=None, alt_loc=None, verbosity=False, merge=False):
1829 """Method for loading structures from a PDB file. 1830 1831 @param file_path: The full path of the PDB file. 1832 @type file_path: str 1833 @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. 1834 @type read_mol: None, int, or list of int 1835 @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. 1836 @type set_mol_name: None, str, or list of str 1837 @keyword read_model: The PDB model to extract from the file. If set to None, then all models will be loaded. 1838 @type read_model: None, int, or list of int 1839 @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. 1840 @type set_model_num: None, int, or list of int 1841 @keyword alt_loc: The PDB ATOM record 'Alternate location indicator' field value to select which coordinates to use. 1842 @type alt_loc: str or None 1843 @keyword verbosity: A flag which if True will cause messages to be printed. 1844 @type verbosity: bool 1845 @keyword merge: A flag which if set to True will try to merge the PDB structure into the currently loaded structures. 1846 @type merge: bool 1847 @return: The status of the loading of the PDB file. 1848 @rtype: bool 1849 """ 1850 1851 # Initial printout. 1852 if verbosity: 1853 print("\nInternal relax PDB parser.") 1854 1855 # Test if the file exists. 1856 if not access(file_path, F_OK): 1857 # Exit indicating failure. 1858 return False 1859 1860 # Separate the file name and path. 1861 path, file = os.path.split(file_path) 1862 1863 # The absolute path. 1864 path_abs = abspath(curdir) + sep + path 1865 1866 # Convert the structure reading args into lists. 1867 if read_mol and not isinstance(read_mol, list): 1868 read_mol = [read_mol] 1869 if set_mol_name and not isinstance(set_mol_name, list): 1870 set_mol_name = [set_mol_name] 1871 if read_model and not isinstance(read_model, list): 1872 read_model = [read_model] 1873 if set_model_num and not isinstance(set_model_num, list): 1874 set_model_num = [set_model_num] 1875 1876 # Open the PDB file. 1877 pdb_file = open_read_file(file_path, verbosity=verbosity) 1878 pdb_lines = pdb_file.readlines() 1879 pdb_file.close() 1880 1881 # Check for empty files. 1882 if pdb_lines == []: 1883 raise RelaxError("The PDB file is empty.") 1884 1885 # Pre-process the lines, fixing PDB violations. 1886 pdb_lines = self._validate_records(pdb_lines) 1887 1888 # Process the different sections. 1889 pdb_lines = self._parse_pdb_title(pdb_lines) 1890 pdb_lines = self._parse_pdb_prim_struct(pdb_lines) 1891 pdb_lines = self._parse_pdb_hetrogen(pdb_lines) 1892 pdb_lines = self._parse_pdb_ss(pdb_lines) 1893 pdb_lines = self._parse_pdb_connectivity_annotation(pdb_lines) 1894 pdb_lines = self._parse_pdb_misc(pdb_lines) 1895 pdb_lines = self._parse_pdb_transform(pdb_lines) 1896 1897 # Loop over all models in the PDB file. 1898 model_index = 0 1899 orig_model_num = [] 1900 mol_conts = [] 1901 for model_num, model_records in self._parse_pdb_coord(pdb_lines): 1902 # Only load the desired model. 1903 if read_model and model_num not in read_model: 1904 continue 1905 1906 # Store the original model number. 1907 orig_model_num.append(model_num) 1908 1909 # Loop over the molecules of the model. 1910 mol_conts.append([]) 1911 mol_index = 0 1912 orig_mol_num = [] 1913 new_mol_name = [] 1914 for mol_num, mol_records in self._parse_mols(model_records): 1915 # Only load the desired model. 1916 if read_mol and mol_num not in read_mol: 1917 continue 1918 1919 # Set the target molecule name. 1920 if set_mol_name: 1921 new_mol_name.append(set_mol_name[mol_index]) 1922 else: 1923 # Number of structures already present for the model. 1924 num_struct = 0 1925 if self.structural_data != None and len(self.structural_data) and (not set_model_num or (model_index <= len(set_model_num) and set_model_num[model_index] == self.structural_data[0].num)): 1926 num_struct = len(self.structural_data[0].mol) 1927 1928 # Set the name to the file name plus the structure number. 1929 new_mol_name.append(file_root(file) + '_mol' + repr(mol_num+num_struct)) 1930 1931 # Store the original mol number. 1932 orig_mol_num.append(mol_num) 1933 1934 # Generate the molecule container. 1935 mol = MolContainer() 1936 1937 # Fill the molecular data object. 1938 mol.fill_object_from_pdb(mol_records, alt_loc_select=alt_loc) 1939 1940 # Store the molecule container. 1941 mol_conts[model_index].append(mol) 1942 1943 # Increment the molecule index. 1944 mol_index = mol_index + 1 1945 1946 # Increment the model index. 1947 model_index = model_index + 1 1948 1949 # No data, so throw a warning and exit. 1950 if not len(mol_conts): 1951 warn(RelaxWarning("No structural data could be read from the file '%s'." % file_path)) 1952 return False 1953 1954 # Create the structural data data structures. 1955 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, file_path_abs=path_abs, merge=merge, verbosity=verbosity) 1956 1957 # Loading worked. 1958 return True
1959 1960
1961 - def load_xyz(self, file_path, read_mol=None, set_mol_name=None, read_model=None, set_model_num=None, verbosity=False):
1962 """Method for loading structures from a XYZ file. 1963 1964 @param file_path: The full path of the XYZ file. 1965 @type file_path: str 1966 @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. 1967 @type read_mol: None, int, or list of int 1968 @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. 1969 @type set_mol_name: None, str, or list of str 1970 @keyword read_model: The XYZ model to extract from the file. If set to None, then all models will be loaded. 1971 @type read_model: None, int, or list of int 1972 @keyword set_model_num: Set the model number of the loaded molecule. If set to None, then the XYZ model numbers will be preserved, if they exist. 1973 @type set_model_num: None, int, or list of int 1974 @keyword verbosity: A flag which if True will cause messages to be printed. 1975 @type verbosity: bool 1976 @return: The status of the loading of the XYZ file. 1977 @rtype: bool 1978 """ 1979 1980 # Initial printout. 1981 if verbosity: 1982 print("\nInternal relax XYZ parser.") 1983 1984 # Test if the file exists. 1985 if not access(file_path, F_OK): 1986 # Exit indicating failure. 1987 return False 1988 1989 # Separate the file name and path. 1990 path, file = os.path.split(file_path) 1991 1992 # The absolute path. 1993 path_abs = abspath(curdir) + sep + path 1994 1995 # Convert the structure reading args into lists. 1996 if read_mol and not isinstance(read_mol, list): 1997 read_mol = [read_mol] 1998 if set_mol_name and not isinstance(set_mol_name, list): 1999 set_mol_name = [set_mol_name] 2000 if read_model and not isinstance(read_model, list): 2001 read_model = [read_model] 2002 if set_model_num and not isinstance(set_model_num, list): 2003 set_model_num = [set_model_num] 2004 2005 # Loop over all models in the XYZ file. 2006 mol_index = 0 2007 model_index = 0 2008 xyz_model_increment = 0 2009 orig_model_num = [] 2010 mol_conts = [] 2011 orig_mol_num = [] 2012 new_mol_name = [] 2013 for model_records in self._parse_models_xyz(file_path, verbosity=verbosity): 2014 # Increment the xyz_model_increment 2015 xyz_model_increment = xyz_model_increment +1 2016 2017 # Only load the desired model. 2018 if read_model and xyz_model_increment not in read_model: 2019 continue 2020 2021 # Store the original model number. 2022 orig_model_num.append(model_index) 2023 2024 # Loop over the molecules of the model. 2025 if read_mol and mol_index not in read_mol: 2026 continue 2027 2028 # Set the target molecule name. 2029 if set_mol_name: 2030 new_mol_name.append(set_mol_name[mol_index]) 2031 else: 2032 if mol_index == 0: 2033 #Set the name to the file name plus the structure number. 2034 new_mol_name.append(file_root(file) + '_mol' + repr(mol_index+1)) 2035 2036 # Store the original mol number. 2037 orig_mol_num.append(mol_index) 2038 2039 # Generate the molecule container. 2040 mol = MolContainer() 2041 2042 # Fill the molecular data object. 2043 mol.fill_object_from_xyz(model_records) 2044 2045 # Store the molecule container. 2046 mol_conts.append([]) 2047 mol_conts[model_index].append(mol) 2048 2049 # Increment the molecule index. 2050 mol_index = mol_index + 1 2051 2052 # Increment the model index. 2053 model_index = model_index + 1 2054 2055 # Create the structural data data structures. 2056 orig_mol_num = [0] 2057 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, file_path_abs=path_abs, verbosity=verbosity) 2058 2059 # Loading worked. 2060 return True
2061 2062
2063 - def model_loop(self, model=None):
2064 """Generator method for looping over the models in numerical order. 2065 2066 @keyword model: Limit the loop to a single number. 2067 @type model: int 2068 @return: The model structural object. 2069 @rtype: ModelContainer container 2070 """ 2071 2072 # A single model. 2073 if model: 2074 for i in range(len(self.structural_data)): 2075 if self.structural_data[i].num == model: 2076 yield self.structural_data[i] 2077 2078 # All models. 2079 else: 2080 # The models. 2081 model_nums = [] 2082 for i in range(len(self.structural_data)): 2083 if self.structural_data[i].num != None: 2084 model_nums.append(self.structural_data[i].num) 2085 2086 # Sort. 2087 if model_nums: 2088 model_nums.sort() 2089 2090 # Loop over the models in order. 2091 for model_num in model_nums: 2092 # Find the model. 2093 for i in range(len(self.structural_data)): 2094 # Yield the model. 2095 if self.structural_data[i].num == model_num: 2096 yield self.structural_data[i] 2097 2098 # No models, so just yield the single container. 2099 if not model_nums: 2100 yield self.structural_data[0]
2101 2102
2103 - def num_models(self):
2104 """Method for returning the number of models. 2105 2106 @return: The number of models in the structural object. 2107 @rtype: int 2108 """ 2109 2110 return len(self.structural_data)
2111 2112
2113 - def num_molecules(self):
2114 """Method for returning the number of molecules. 2115 2116 @return: The number of molecules in the structural object. 2117 @rtype: int 2118 """ 2119 2120 # No data. 2121 if self.empty(): 2122 return 0 2123 2124 # Validate the structural object. 2125 self.validate() 2126 2127 # Return the number. 2128 return len(self.structural_data[0].mol)
2129 2130
2131 - def pack_structs(self, data_matrix, orig_model_num=None, set_model_num=None, orig_mol_num=None, set_mol_name=None, file_name=None, file_path=None, file_path_abs=None, verbosity=1, merge=False):
2132 """From the given structural data, expand the structural data data structure. 2133 2134 @param data_matrix: A matrix of structural objects. 2135 @type data_matrix: list of lists of structural objects 2136 @keyword orig_model_num: The original model numbers (for storage). 2137 @type orig_model_num: list of int 2138 @keyword set_model_num: The new model numbers (for model renumbering). 2139 @type set_model_num: list of int 2140 @keyword orig_mol_num: The original molecule numbers (for storage). 2141 @type orig_mol_num: list of int 2142 @keyword set_mol_name: The new molecule names. 2143 @type set_mol_name: list of str 2144 @keyword file_name: The name of the file from which the molecular data has been extracted. 2145 @type file_name: None or str 2146 @keyword file_path: The full path to the file specified by 'file_name'. 2147 @type file_path: None or str 2148 @keyword file_path_abs: The absolute path to the file specified by 'file_name'. This is a fallback mechanism in case results or save files are located somewhere other than the working directory. 2149 @type file_path_abs: None or str 2150 @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. 2151 @type verbosity: int 2152 @keyword merge: A flag which if set to True will try to merge the structure into the currently loaded structures. 2153 @type merge: bool 2154 """ 2155 2156 # Test the number of models. 2157 if len(orig_model_num) != len(data_matrix): 2158 raise RelaxError("Structural data mismatch, %s original models verses %s in the structural data." % (len(orig_model_num), len(data_matrix))) 2159 2160 # Test the number of molecules. 2161 if len(orig_mol_num) != len(data_matrix[0]): 2162 raise RelaxError("Structural data mismatch, %s original molecules verses %s in the structural data." % (len(orig_mol_num), len(data_matrix[0]))) 2163 2164 # Model numbers do not change. 2165 if not set_model_num: 2166 set_model_num = orig_model_num 2167 2168 # Test the model mapping. 2169 if len(set_model_num) != len(data_matrix): 2170 raise RelaxError("Failure of the mapping of new model numbers, %s new model numbers verses %s models in the structural data." % (len(set_model_num), len(data_matrix))) 2171 2172 # Test the molecule mapping. 2173 if len(set_mol_name) != len(data_matrix[0]): 2174 raise RelaxError("Failure of the mapping of new molecule names, %s new molecule names verses %s molecules in the structural data." % (len(set_mol_name), len(data_matrix[0]))) 2175 2176 # Test that the target models and structures are absent, and get the already existing model numbers. 2177 for i in range(len(set_model_num)): 2178 # Merging flag is set, so skip the checks. 2179 if merge: 2180 continue 2181 2182 # A new model, so no need to check. 2183 if not set_model_num[i] in self.structural_data.current_models: 2184 continue 2185 2186 # Loop over the structures. 2187 index = self.structural_data.current_models.index(set_model_num[i]) 2188 for j in range(len(self.structural_data[index].mol)): 2189 if self.structural_data[index].mol[j].mol_name in set_mol_name: 2190 raise RelaxError("The molecule '%s' of model %s already exists." % (self.structural_data[i].mol[j].mol_name, self.structural_data[i].num)) 2191 2192 # Loop over the models. 2193 for i in range(len(set_model_num)): 2194 # The model doesn't currently exist. 2195 if set_model_num[i] not in self.structural_data.current_models: 2196 # Create the model. 2197 self.structural_data.add_item(set_model_num[i]) 2198 2199 # Get the model. 2200 model = self.structural_data[-1] 2201 2202 # Otherwise get the pre-existing model. 2203 else: 2204 model = self.structural_data[self.structural_data.current_models.index(set_model_num[i])] 2205 2206 # Loop over the molecules. 2207 for j in range(len(set_mol_name)): 2208 # Printout. 2209 if verbosity: 2210 if merge: 2211 print("Merging with model %s of molecule '%s' (from the original molecule number %s of model %s)" % (set_model_num[i], set_mol_name[j], orig_mol_num[j], orig_model_num[i])) 2212 else: 2213 print("Adding molecule '%s' to model %s (from the original molecule number %s of model %s)" % (set_mol_name[j], set_model_num[i], orig_mol_num[j], orig_model_num[i])) 2214 2215 # The index of the new molecule to add or merge. 2216 index = len(model.mol) 2217 if merge: 2218 index -= 1 2219 2220 # Consistency check. 2221 if model.num != self.structural_data[0].num and self.structural_data[0].mol[index].mol_name != set_mol_name[j]: 2222 raise RelaxError("The new molecule name of '%s' in model %s does not match the corresponding molecule's name of '%s' in model %s." % (set_mol_name[j], set_model_num[i], self.structural_data[0].mol[index].mol_name, self.structural_data[0].num)) 2223 2224 # Pack the structures. 2225 if merge: 2226 mol = model.mol.merge_item(mol_name=set_mol_name[j], mol_cont=data_matrix[i][j]) 2227 else: 2228 mol = model.mol.add_item(mol_name=set_mol_name[j], mol_cont=data_matrix[i][j]) 2229 2230 # Set the molecule name and store the structure file info. 2231 mol.mol_name = set_mol_name[j] 2232 mol.file_name = file_name 2233 mol.file_path = file_path 2234 mol.file_path_abs = file_path_abs 2235 mol.file_mol_num = orig_mol_num[j] 2236 mol.file_model = orig_model_num[i]
2237 2238
2239 - def rotate(self, R=None, origin=None, model=None, atom_id=None):
2240 """Rotate the structural information about the given origin. 2241 2242 @keyword R: The forwards rotation matrix. 2243 @type R: numpy 3D, rank-2 array 2244 @keyword origin: The origin of the rotation. 2245 @type origin: numpy 3D, rank-1 array 2246 @keyword model: The model to rotate. If None, all models will be rotated. 2247 @type model: int 2248 @keyword atom_id: The molecule, residue, and atom identifier string. Only atoms matching this selection will be used. 2249 @type atom_id: str or None 2250 """ 2251 2252 # Generate the selection object. 2253 sel_obj = None 2254 if atom_id: 2255 sel_obj = Selection(atom_id) 2256 2257 # Loop over the models. 2258 for model_cont in self.model_loop(model): 2259 # Loop over the molecules. 2260 for mol in model_cont.mol: 2261 # Skip non-matching molecules. 2262 if sel_obj and not sel_obj.contains_mol(mol.mol_name): 2263 continue 2264 2265 # Loop over the atoms. 2266 for i in range(len(mol.atom_num)): 2267 # Skip non-matching atoms. 2268 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): 2269 continue 2270 2271 # The origin to atom vector. 2272 vect = array([mol.x[i], mol.y[i], mol.z[i]], float64) - origin 2273 2274 # Rotation. 2275 rot_vect = dot(R, vect) 2276 2277 # The new position. 2278 pos = rot_vect + origin 2279 mol.x[i] = pos[0] 2280 mol.y[i] = pos[1] 2281 mol.z[i] = pos[2]
2282 2283
2284 - def target_mol_name(self, set=None, target=None, index=None, mol_num=None, file=None):
2285 """Add the new molecule name to the target data structure. 2286 2287 @keyword set: The list of new molecule names. If not supplied, the names will be generated from the file name. 2288 @type set: None or list of str 2289 @keyword target: The target molecule name data structure to which the new name will be appended. 2290 @type target: list 2291 @keyword index: The molecule index, matching the set argument. 2292 @type index: int 2293 @keyword mol_num: The molecule number. 2294 @type mol_num: int 2295 @keyword file: The name of the file, excluding all directories. 2296 @type file: str 2297 """ 2298 2299 # Set the target molecule name. 2300 if set: 2301 target.append(set[index]) 2302 else: 2303 # Set the name to the file name plus the structure number. 2304 target.append(file_root(file) + '_mol' + repr(mol_num))
2305 2306
2307 - def translate(self, T=None, model=None, atom_id=None):
2308 """Displace the structural information by the given translation vector. 2309 2310 @keyword T: The translation vector. 2311 @type T: numpy 3D, rank-1 array 2312 @keyword model: The model to rotate. If None, all models will be rotated. 2313 @type model: int 2314 @keyword atom_id: The molecule, residue, and atom identifier string. Only atoms matching this selection will be used. 2315 @type atom_id: str or None 2316 """ 2317 2318 # Generate the selection object. 2319 sel_obj = None 2320 if atom_id: 2321 sel_obj = Selection(atom_id) 2322 2323 # Loop over the models. 2324 for model_cont in self.model_loop(model): 2325 # Loop over the molecules. 2326 for mol in model_cont.mol: 2327 # Skip non-matching molecules. 2328 if sel_obj and not sel_obj.contains_mol(mol.mol_name): 2329 continue 2330 2331 # Loop over the atoms. 2332 for i in range(len(mol.atom_num)): 2333 # Skip non-matching atoms. 2334 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): 2335 continue 2336 2337 # Translate. 2338 mol.x[i] = mol.x[i] + T[0] 2339 mol.y[i] = mol.y[i] + T[1] 2340 mol.z[i] = mol.z[i] + T[2]
2341 2342
2343 - def to_xml(self, doc, element):
2344 """Prototype method for converting the structural object to an XML representation. 2345 2346 @param doc: The XML document object. 2347 @type doc: xml.dom.minidom.Document instance 2348 @param element: The element to add the alignment tensors XML element to. 2349 @type element: XML element object 2350 """ 2351 2352 # Create the structural element and add it to the higher level element. 2353 str_element = doc.createElement('structure') 2354 element.appendChild(str_element) 2355 2356 # Set the structural attributes. 2357 str_element.setAttribute('desc', 'Structural information') 2358 2359 # No contents to store, so pack up the structural containers. 2360 if not self.structural_data.is_empty(): 2361 self.structural_data.to_xml(doc, str_element) 2362 2363 # The structural metadata. 2364 metadata = ['helices', 'sheets'] 2365 for name in metadata: 2366 # The metadata does not exist. 2367 if not hasattr(self, name): 2368 continue 2369 2370 # Get the object. 2371 obj = getattr(self, name) 2372 2373 # Create a new element for this object, and add it to the main element. 2374 sub_elem = doc.createElement(name) 2375 str_element.appendChild(sub_elem) 2376 2377 # Add the value to the sub element. 2378 object_to_xml(doc, sub_elem, value=obj) 2379 2380 # The displacement structure. 2381 if hasattr(self, 'displacements'): 2382 # Create an XML element. 2383 disp_element = doc.createElement('displacements') 2384 str_element.appendChild(disp_element) 2385 2386 # Set the attributes. 2387 disp_element.setAttribute('desc', 'The rotational and translational displacements between models') 2388 2389 # Add the displacement data. 2390 self.displacements.to_xml(doc, disp_element)
2391 2392
2393 - def validate_models(self):
2394 """Check that the models are consistent with each other. 2395 2396 This checks that the primary structure is identical between the models. 2397 """ 2398 2399 # Print out. 2400 print("Validating models:") 2401 2402 # Loop over the models. 2403 for i in range(len(self.structural_data)): 2404 # Check the molecules. 2405 if len(self.structural_data[0].mol) != len(self.structural_data[i].mol): 2406 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))) 2407 2408 # Loop over the molecules. 2409 for j in range(len(self.structural_data[i].mol)): 2410 # Alias the molecules. 2411 mol = self.structural_data[i].mol[j] 2412 mol_ref = self.structural_data[0].mol[j] 2413 2414 # Check the names. 2415 if mol.mol_name != mol_ref.mol_name: 2416 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)) 2417 2418 # Loop over the atoms. 2419 for k in range(len(mol.atom_name)): 2420 # Create pseudo-pdb formatted records (with no atomic coordinates). 2421 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]), '') 2422 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]), '') 2423 2424 # Check the atom info. 2425 if atom != atom_ref: 2426 print(atom) 2427 print(atom_ref) 2428 raise RelaxError("The atoms of model %i do not match the first model." % self.structural_data[i].num) 2429 2430 # Final printout. 2431 print("\tAll models are consistent")
2432 2433
2434 - def write_pdb(self, file, model_num=None):
2435 """Method for the creation of a PDB file from the structural data. 2436 2437 A number of PDB records including HET, HETNAM, FORMUL, HELIX, SHEET, HETATM, TER, CONECT, MASTER, and END are created. To create the non-standard residue records HET, HETNAM, and FORMUL, the data structure 'het_data' is created. It is an array of arrays where the first dimension corresponds to a different residue and the second dimension has the elements: 2438 2439 0. Residue number. 2440 1. Residue name. 2441 2. Chain ID. 2442 3. Total number of atoms in the residue. 2443 4. Number of H atoms in the residue. 2444 5. Number of C atoms in the residue. 2445 2446 2447 @param file: The PDB file object. This object must be writable. 2448 @type file: file object 2449 @keyword model_num: The model to place into the PDB file. If not supplied, then all models will be placed into the file. 2450 @type model_num: None or int 2451 """ 2452 2453 # Validate the structural data. 2454 self.validate() 2455 2456 # Print out. 2457 print("\nCreating the PDB records\n") 2458 2459 # Write some initial remarks. 2460 print("REMARK") 2461 pdb_write.remark(file, num=4, remark="This file complies with format v. 3.30, Jul-2011.") 2462 pdb_write.remark(file, num=40, remark=None) 2463 pdb_write.remark(file, num=40, remark="Created using relax (http://www.nmr-relax.com).") 2464 pdb_write.remark(file, num=40, remark=None) 2465 if RELAX_VERSION: 2466 pdb_write.remark(file, num=40, remark="relax version %s." % RELAX_VERSION) 2467 pdb_write.remark(file, num=40, remark="Created on %s." % asctime()) 2468 num_remark = 2 2469 2470 # Determine if model records will be created. 2471 model_records = False 2472 for model in self.model_loop(): 2473 if hasattr(model, 'num') and model.num != None: 2474 model_records = True 2475 2476 2477 #################### 2478 # Hetrogen section # 2479 #################### 2480 2481 # Initialise the hetrogen info array. 2482 het_data = [] 2483 het_data_coll = [] 2484 2485 # Loop over the molecules of the first model. 2486 index = 0 2487 for mol in self.structural_data[0].mol: 2488 # Check the validity of the data. 2489 self._validate_data_arrays(mol) 2490 2491 # Append an empty array for this molecule. 2492 het_data.append([]) 2493 2494 # Collect the non-standard residue info. 2495 for i in range(len(mol.atom_name)): 2496 # Skip non-HETATM records and HETATM records with no residue info. 2497 if mol.pdb_record[i] != 'HETATM' or mol.res_name[i] == None: 2498 continue 2499 2500 # If the residue is not already stored initialise a new het_data element. 2501 # (residue number, residue name, chain ID, number of atoms, atom count array). 2502 if not het_data[index] or not mol.res_num[i] == het_data[index][-1][0]: 2503 het_data[index].append([mol.res_num[i], mol.res_name[i], CHAIN_ID_LIST[index], 0, []]) 2504 2505 # Catch missing chain_ids. 2506 if het_data[index][-1][2] == None: 2507 het_data[index][-1][2] = '' 2508 2509 # Total atom count. 2510 het_data[index][-1][3] = het_data[index][-1][3] + 1 2511 2512 # Find if the atom has already a count entry. 2513 entry = False 2514 for j in range(len(het_data[index][-1][4])): 2515 if mol.element[i] == het_data[index][-1][4][j][0]: 2516 entry = True 2517 2518 # Create a new specific atom count entry. 2519 if not entry: 2520 het_data[index][-1][4].append([mol.element[i], 0]) 2521 2522 # Increment the specific atom count. 2523 for j in range(len(het_data[index][-1][4])): 2524 if mol.element[i] == het_data[index][-1][4][j][0]: 2525 het_data[index][-1][4][j][1] = het_data[index][-1][4][j][1] + 1 2526 2527 # Create the collective hetrogen info data structure. 2528 for i in range(len(het_data[index])): 2529 # Find the entry in the collective structure. 2530 found = False 2531 for j in range(len(het_data_coll)): 2532 # Different residue numbers. 2533 if het_data[index][i][0] != het_data_coll[j][0]: 2534 continue 2535 2536 # Different chain IDs. 2537 if het_data_coll[j][2] != het_data[index][i][2]: 2538 continue 2539 2540 # Change the flag. 2541 found = True 2542 2543 # If there is no match, add the new residue to the collective. 2544 if not found: 2545 het_data_coll.append(het_data[index][i]) 2546 2547 # Increment the molecule index. 2548 index = index + 1 2549 2550 2551 # The HET records. 2552 ################## 2553 2554 # Print out. 2555 print("HET") 2556 2557 # Write the HET records. 2558 for het in het_data_coll: 2559 pdb_write.het(file, het_id=het[1], chain_id=het[2], seq_num=het[0], num_het_atoms=het[3]) 2560 2561 2562 # The HETNAM records. 2563 ##################### 2564 2565 # Print out. 2566 print("HETNAM") 2567 2568 # Loop over the non-standard residues. 2569 residues = [] 2570 for het in het_data_coll: 2571 # Test if the residue HETNAM record as already been written (otherwise store its name). 2572 if het[1] in residues: 2573 continue 2574 else: 2575 residues.append(het[1]) 2576 2577 # Get the chemical name. 2578 chemical_name = self._get_chemical_name(het[1]) 2579 if not chemical_name: 2580 chemical_name = 'Unknown' 2581 2582 # Write the HETNAM records. 2583 pdb_write.hetnam(file, het_id=het[1], text=chemical_name) 2584 2585 2586 # The FORMUL records. 2587 ##################### 2588 2589 # Print out. 2590 print("FORMUL") 2591 2592 # Loop over the non-standard residues and generate and write the chemical formula. 2593 residues = [] 2594 for i in range(len(het_data_coll)): 2595 # Alias. 2596 het = het_data_coll[i] 2597 2598 # Test if the residue HETNAM record as already been written (otherwise store its name). 2599 if het[1] in residues: 2600 continue 2601 else: 2602 residues.append(het[1]) 2603 2604 # Initialise the chemical formula. 2605 formula = '' 2606 2607 # Loop over the atoms. 2608 for atom_count in het[4]: 2609 formula = formula + atom_count[0] + repr(atom_count[1]) 2610 2611 # The FORMUL record (chemical formula). 2612 pdb_write.formul(file, comp_num=i+1, het_id=het[1], text=formula) 2613 2614 2615 ############################### 2616 # Secondary structure section # 2617 ############################### 2618 2619 # The HELIX records. 2620 #################### 2621 2622 if hasattr(self, 'helices') and len(self.helices): 2623 # Printout. 2624 print("HELIX") 2625 2626 # Loop over and unpack the helix data. 2627 index = 1 2628 for helix_id, init_chain_id, init_res_name, init_seq_num, end_chain_id, end_res_name, end_seq_num, helix_class, length in self.helices: 2629 pdb_write.helix(file, ser_num=index, helix_id=helix_id, init_chain_id=init_chain_id, init_res_name=init_res_name, init_seq_num=init_seq_num, end_chain_id=end_chain_id, end_res_name=end_res_name, end_seq_num=end_seq_num, helix_class=helix_class, length=length) 2630 index += 1 2631 2632 # The SHEET records. 2633 #################### 2634 2635 if hasattr(self, 'sheets') and len(self.sheets): 2636 # Printout. 2637 print("SHEET") 2638 2639 # Loop over and unpack the helix data. 2640 index = 1 2641 for strand, sheet_id, num_strands, init_res_name, init_chain_id, init_seq_num, init_icode, end_res_name, end_chain_id, end_seq_num, end_icode, sense, cur_atom, cur_res_name, cur_chain_id, cur_res_seq, cur_icode, prev_atom, prev_res_name, prev_chain_id, prev_res_seq, prev_icode in self.sheets: 2642 pdb_write.sheet(file, strand=strand, sheet_id=sheet_id, num_strands=num_strands, init_res_name=init_res_name, init_chain_id=init_chain_id, init_seq_num=init_seq_num, init_icode=init_icode, end_res_name=end_res_name, end_chain_id=end_chain_id, end_seq_num=end_seq_num, end_icode=end_icode, sense=sense, cur_atom=cur_atom, cur_res_name=cur_res_name, cur_chain_id=cur_chain_id, cur_res_seq=cur_res_seq, cur_icode=cur_icode, prev_atom=prev_atom, prev_res_name=prev_res_name, prev_chain_id=prev_chain_id, prev_res_seq=prev_res_seq, prev_icode=prev_icode) 2643 index += 1 2644 2645 2646 ###################### 2647 # Coordinate section # 2648 ###################### 2649 2650 # Loop over the models. 2651 for model in self.model_loop(model_num): 2652 # Initialise record counts. 2653 num_hetatm = 0 2654 num_atom = 0 2655 num_ter = 0 2656 ser_num = 1 2657 2658 # MODEL record, for multiple models. 2659 #################################### 2660 2661 if model_records: 2662 # Print out. 2663 print("\nMODEL %s" % model.num) 2664 2665 # Write the model record. 2666 pdb_write.model(file, serial=model.num) 2667 2668 2669 # Add the atomic coordinate records (ATOM, HETATM, and TER). 2670 ############################################################ 2671 2672 # Loop over the molecules. 2673 index = 0 2674 atom_serial = 0 2675 for mol in model.mol: 2676 # Print out. 2677 print("ATOM, HETATM, TER") 2678 2679 # Loop over the atomic data. 2680 atom_record = False 2681 for i in range(len(mol.atom_name)): 2682 # Write the ATOM record. 2683 if mol.pdb_record[i] in [None, 'ATOM']: 2684 atom_record = True 2685 2686 # The atom number. 2687 atom_serial += 1 2688 2689 # Handle the funky atom name alignment. From the PDB format documents: 2690 # "Alignment of one-letter atom name such as C starts at column 14, while two-letter atom name such as FE starts at column 13." 2691 if len(mol.atom_name[i]) == 1: 2692 atom_name = " %s" % mol.atom_name[i] 2693 else: 2694 atom_name = "%s" % mol.atom_name[i] 2695 2696 # Write out. 2697 pdb_write.atom(file, serial=atom_serial, name=atom_name, res_name=mol.res_name[i], chain_id=CHAIN_ID_LIST[index], res_seq=mol.res_num[i], x=mol.x[i], y=mol.y[i], z=mol.z[i], occupancy=1.0, temp_factor=0, element=mol.element[i]) 2698 num_atom += 1 2699 ser_num += 1 2700 2701 # Info for the TER record. 2702 ter_num = atom_serial + 1 2703 ter_name = mol.res_name[i] 2704 ter_chain_id = CHAIN_ID_LIST[index] 2705 ter_res_num = mol.res_num[i] 2706 2707 # Finish the ATOM section with the TER record. 2708 if atom_record: 2709 pdb_write.ter(file, serial=ser_num, res_name=ter_name, chain_id=CHAIN_ID_LIST[index], res_seq=ter_res_num) 2710 num_ter += 1 2711 ser_num += 1 2712 2713 # Loop over the atomic data. 2714 count_shift = False 2715 for i in range(len(mol.atom_name)): 2716 # Write the HETATM record. 2717 if mol.pdb_record[i] == 'HETATM': 2718 # The atom number. 2719 atom_serial += 1 2720 2721 # Increment the atom number if a TER record was created. 2722 if atom_record and atom_serial == ter_num: 2723 count_shift = True 2724 if atom_record and count_shift: 2725 atom_serial += 1 2726 2727 # Write out. 2728 pdb_write.hetatm(file, serial=ser_num, name=self._translate(mol.atom_name[i]), res_name=mol.res_name[i], chain_id=CHAIN_ID_LIST[index], res_seq=mol.res_num[i], x=mol.x[i], y=mol.y[i], z=mol.z[i], occupancy=1.0, temp_factor=0.0, element=mol.element[i]) 2729 num_hetatm += 1 2730 ser_num += 1 2731 2732 # Increment the molecule index. 2733 index += 1 2734 2735 2736 # ENDMDL record, for multiple structures. 2737 ######################################## 2738 2739 if model_records: 2740 print("ENDMDL") 2741 pdb_write.endmdl(file) 2742 2743 2744 # Create the CONECT records. 2745 ############################ 2746 2747 # Print out. 2748 print("CONECT") 2749 2750 # Initialise record counts. 2751 num_conect = 0 2752 2753 # The per molecule incremented atom counts. 2754 atom_counts = [0] 2755 index = 0 2756 for mol in self.structural_data[0].mol: 2757 if index == 0: 2758 atom_counts.append(len(mol.atom_name)) 2759 else: 2760 atom_counts.append(atom_counts[index] + len(mol.atom_name)) 2761 index += 1 2762 2763 # Loop over the molecules of the first model. 2764 index = 0 2765 for mol in self.structural_data[0].mol: 2766 # Loop over the atoms. 2767 for i in range(len(mol.atom_name)): 2768 # No bonded atoms, hence no CONECT record is required. 2769 if not len(mol.bonded[i]): 2770 continue 2771 2772 # Initialise some data structures. 2773 flush = 0 2774 bonded_index = 0 2775 bonded = ['', '', '', ''] 2776 bonded_shifted = ['', '', '', ''] 2777 2778 # Loop over the bonded atoms. 2779 for j in range(len(mol.bonded[i])): 2780 # End of the array, hence create the CONECT record in this iteration. 2781 if j == len(mol.bonded[i])-1: 2782 flush = True 2783 2784 # Only four covalently bonded atoms allowed in one CONECT record. 2785 if bonded_index == 3: 2786 flush = True 2787 2788 # Get the bonded atom index. 2789 bonded[bonded_index] = mol.bonded[i][j] 2790 2791 # Increment the bonded_index value. 2792 bonded_index = bonded_index + 1 2793 2794 # Generate the CONECT record and increment the counter. 2795 if flush: 2796 # Convert the atom indices to atom numbers. 2797 for k in range(4): 2798 if bonded[k] != '': 2799 if mol.atom_num[bonded[k]] != None: 2800 bonded[k] = mol.atom_num[bonded[k]] 2801 else: 2802 bonded[k] = bonded[k] + 1 2803 bonded_shifted[k] = bonded[k] + atom_counts[index] 2804 2805 # Write the CONECT record. 2806 pdb_write.conect(file, serial=i+1+atom_counts[index], bonded1=bonded_shifted[0], bonded2=bonded_shifted[1], bonded3=bonded_shifted[2], bonded4=bonded_shifted[3]) 2807 2808 # Reset the flush flag, the bonded atom count, and the bonded atom names. 2809 flush = False 2810 bonded_index = 0 2811 bonded = ['', '', '', ''] 2812 2813 # Increment the CONECT record count. 2814 num_conect = num_conect + 1 2815 2816 # Increment the molecule index. 2817 index += 1 2818 2819 2820 # MASTER record. 2821 ################ 2822 2823 print("\nMASTER") 2824 pdb_write.master(file, num_het=len(het_data_coll), num_coord=num_atom+num_hetatm, num_ter=num_ter, num_conect=num_conect) 2825 2826 2827 # END. 2828 ###### 2829 2830 print("END") 2831 pdb_write.end(file)
2832 2833
2834 - def validate(self):
2835 """Check the integrity of the structural data. 2836 2837 The number of molecules must be the same in all models. 2838 """ 2839 2840 # Reference number of molecules. 2841 num_mols = len(self.structural_data[0].mol) 2842 2843 # Loop over all other models. 2844 for i in range(1, len(self.structural_data)): 2845 # Model alias. 2846 model_i = self.structural_data[i] 2847 2848 # Size check. 2849 if num_mols != len(model_i.mol): 2850 raise RelaxError("The structural object is not valid - the number of molecules is not the same for all models.") 2851 2852 # Molecule name check. 2853 for j in range(len(model_i.mol)): 2854 if model_i.mol[j].mol_name != self.structural_data[0].mol[j].mol_name: 2855 raise RelaxError("The molecule name '%s' of model %s does not match the corresponding molecule '%s' of model %s." % (model_i.mol[j].mol_name, model_i.num, self.structural_data[0].mol[j].mol_name, self.structural_data[0].num))
2856