lib.structure.internal.object

62 """The internal relax structural data object. 63 64 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. 65 """ 66

67 - def __init__(self):

68 """Initialise the structural object.""" 69 70 # Initialise the empty model list. 71 self.structural_data = ModelList()

72 73

74 - def _bonded_atom(self, attached_atom, index, mol):

75 """Find the atom named attached_atom directly bonded to the atom located at the index. 76 77 @param attached_atom: The name of the attached atom to return. 78 @type attached_atom: str 79 @param index: The index of the atom which the attached atom is attached to. 80 @type index: int 81 @param mol: The molecule container. 82 @type mol: MolContainer instance 83 @return: A tuple of information about the bonded atom. 84 @rtype: tuple consisting of the atom number (int), atom name (str), element name (str), and atomic position (Numeric array of len 3) 85 """ 86 87 # Init. 88 bonded_found = False 89 90 # No bonded atoms, so determine the connectivities. 91 if not mol.bonded[index]: 92 # Determine the molecule type if needed. 93 if not hasattr(mol, 'type'): 94 self._mol_type(mol) 95 96 # Protein. 97 if mol.type == 'protein': 98 self._protein_connect(mol) 99 100 # Find everything within 2 Angstroms and say they are bonded. 101 else: 102 self._find_bonded_atoms(index, mol, radius=2) 103 104 # Loop over the bonded atoms. 105 matching_list = [] 106 for bonded_index in mol.bonded[index]: 107 if regex.search(mol.atom_name[bonded_index], attached_atom): 108 matching_list.append(bonded_index) 109 num_attached = len(matching_list) 110 111 # Problem. 112 if num_attached > 1: 113 # Get the atom names. 114 matching_names = [] 115 for i in matching_list: 116 matching_names.append(mol.atom_name[i]) 117 118 # Return nothing but a warning. 119 return None, None, None, None, None, 'More than one attached atom found: ' + repr(matching_names) 120 121 # No attached atoms. 122 if num_attached == 0: 123 return None, None, None, None, None, "No attached atom could be found" 124 125 # The bonded atom info. 126 index = matching_list[0] 127 bonded_num = mol.atom_num[index] 128 bonded_name = mol.atom_name[index] 129 element = mol.element[index] 130 pos = [mol.x[index], mol.y[index], mol.z[index]] 131 attached_name = mol.atom_name[index] 132 133 # Return the information. 134 return bonded_num, bonded_name, element, pos, attached_name, None

135 136

137 - def _find_bonded_atoms(self, index, mol, radius=1.2):

138 """Find all atoms within a sphere and say that they are attached to the central atom. 139 140 The found atoms will be added to the 'bonded' data structure. 141 142 143 @param index: The index of the central atom. 144 @type index: int 145 @param mol: The molecule container. 146 @type mol: MolContainer instance 147 """ 148 149 # Central atom info. 150 centre = array([mol.x[index], mol.y[index], mol.z[index]], float64) 151 152 # Atom loop. 153 dist_list = [] 154 connect_list = {} 155 element_list = {} 156 for i in range(len(mol.atom_num)): 157 # Skip proton to proton bonds! 158 if mol.element[index] == 'H' and mol.element[i] == 'H': 159 continue 160 161 # The atom's position. 162 pos = array([mol.x[i], mol.y[i], mol.z[i]], float64) 163 164 # The distance from the centre. 165 dist = linalg.norm(centre-pos) 166 167 # The atom is within the radius. 168 if dist < radius: 169 # Store the distance. 170 dist_list.append(dist) 171 172 # Store the atom index. 173 connect_list[dist] = i 174 175 # Store the element type. 176 element_list[dist] = mol.element[i] 177 178 # The maximum number of allowed covalent bonds. 179 max_conn = 1000 # Ridiculous default! 180 if mol.element[index] == 'H': 181 max_conn = 1 182 elif mol.element[index] == 'O': 183 max_conn = 2 184 elif mol.element[index] == 'N': 185 max_conn = 3 186 elif mol.element[index] == 'C': 187 max_conn = 4 188 189 # Sort. 190 dist_list.sort() 191 192 # Loop over the max number of connections (or the number of connected atoms, if less). 193 for i in range(min(max_conn, len(dist_list))): 194 mol.atom_connect(index, connect_list[dist_list[i]])

195 196

197 - def _get_chemical_name(self, hetID):

198 """Return the chemical name corresponding to the given residue ID. 199 200 The following names are currently returned:: 201 ___________________________________________________________ 202 | | | 203 | hetID | Chemical name | 204 |________|________________________________________________| 205 | | | 206 | TNS | Tensor | 207 | COM | Centre of mass | 208 | AXS | Tensor axes | 209 | SIM | Monte Carlo simulation tensor axes | 210 | PIV | Pivot point | 211 | AVE | Average vector | 212 | RTX | Axis of the rotor geometric object | 213 | RTB | Propeller blades of the rotor geometric object | 214 | RTL | Labels for the rotor geometric object | 215 | CON | Cone object | 216 | CNC | Apex or centre of the cone geometric object | 217 | CNX | Axis of the cone geometric object | 218 | CNE | Edge of the cone geometric object | 219 | AXE | The axis geometric object | 220 | TLE | The title for the object | 221 |________|________________________________________________| 222 223 For any other residues, no description is returned. 224 225 @param hetID: The residue ID. 226 @type hetID: str 227 @return: The chemical name. 228 @rtype: str or None 229 """ 230 231 # Translation table. 232 table = { 233 "TNS": "Tensor", 234 "COM": "Centre of mass", 235 "AXS": "Tensor axes", 236 "SIM": "Monte Carlo simulation tensor axes", 237 "PIV": "Pivot point", 238 "AVE": "Average vector", 239 "RTX": "Axis of the rotor geometric object", 240 "RTB": "Propeller blades of the rotor geometric object", 241 "RTL": "Labels for the rotor geometric object", 242 "CON": "Cone object", 243 "CNC": "Apex or centre of the cone geometric object", 244 "CNX": "Axis of the cone geometric object", 245 "CNE": "Edge of the cone geometric object", 246 "AXE": "The axis geometric object", 247 "TLE": "The title for the object", 248 } 249 250 # Return the description, if one exists. 251 if hetID in table: 252 return table[hetID]

253 254

255 - def _parse_models_gaussian(self, file_path, verbosity=1):

256 """Generator function for looping over the models in the Gaussian log file. 257 258 @param file_path: The full path of the Gaussian log file. 259 @type file_path: str 260 @return: The model number and all the records for that model. 261 @rtype: tuple of int and array of str 262 @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. 263 @type verbosity: int 264 """ 265 266 # Open the file. 267 file = open_read_file(file_path, verbosity=verbosity) 268 lines = file.readlines() 269 file.close() 270 271 # Check for empty files. 272 if lines == []: 273 raise RelaxError("The Gaussian log file is empty.") 274 275 # Init. 276 found = False 277 str_index = 0 278 total_atom = 0 279 model = 0 280 records = [] 281 282 # Loop over the data. 283 for i in range(len(lines)): 284 # Found a structure. 285 if search("Standard orientation", lines[i]): 286 found = True 287 str_index = 0 288 continue 289 290 # End of the model. 291 if found and str_index > 4 and search("---------", lines[i]): 292 # Yield the info 293 yield records 294 295 # Reset. 296 records = [] 297 found = False 298 299 # Not a structure line. 300 if not found: 301 continue 302 303 # Append the line as a record of the model. 304 records.append(lines[i]) 305 306 # Increment the structure line index. 307 str_index += 1

308 309

310 - def _parse_pdb_connectivity_annotation(self, lines):

311 """Loop over and parse the PDB connectivity annotation records. 312 313 These are the records identified in the U{PDB version 3.30 documentation<http://www.wwpdb.org/documentation/file-format/format33/sect6.html>}. 314 315 316 @param lines: The lines of the PDB file excluding the sections prior to the connectivity annotation section. 317 @type lines: list of str 318 @return: The remaining PDB lines with the connectivity annotation records stripped. 319 @rtype: list of str 320 """ 321 322 # The ordered list of record names in the connectivity annotation section. 323 records = [ 324 'SSBOND', 325 'LINK ', 326 'CISPEP' 327 ] 328 329 # Loop over the lines. 330 for i in range(len(lines)): 331 # No match, therefore assume to be out of the connectivity annotation section. 332 if lines[i][:6] not in records: 333 break 334 335 # Return the remaining lines. 336 return lines[i:]

337 338

339 - def _parse_pdb_coord(self, lines):

340 """Generator function for looping over the models in the PDB file. 341 342 These are the records identified in the PDB version 3.30 documentation at U{http://www.wwpdb.org/documentation/file-format/format33/sect9.html}. 343 344 345 @param lines: The lines of the coordinate section. 346 @type lines: list of str 347 @return: The model number and all the records for that model. 348 @rtype: tuple of int and array of str 349 """ 350 351 # Init. 352 model = None 353 records = [] 354 355 # Loop over the data. 356 for i in range(len(lines)): 357 # A new model record. 358 if lines[i][:5] == 'MODEL': 359 try: 360 model = int(lines[i].split()[1]) 361 except: 362 raise RelaxError("The MODEL record " + repr(lines[i]) + " is corrupt, cannot read the PDB file.") 363 364 # Skip all records prior to the first ATOM or HETATM record. 365 if not (lines[i][:4] == 'ATOM' or lines[i][:6] == 'HETATM') and not len(records): 366 continue 367 368 # End of the model. 369 if lines[i][:6] == 'ENDMDL': 370 # Yield the info. 371 yield model, records 372 373 # Reset the records. 374 records = [] 375 376 # Skip the rest of this loop. 377 continue 378 379 # Append the line as a record of the model. 380 records.append(lines[i]) 381 382 # If records is not empty then there are no models, so yield the lot. 383 if len(records): 384 yield model, records

385 386

387 - def _parse_pdb_hetrogen(self, lines):

388 """Loop over and parse the PDB hetrogen records. 389 390 These are the records identified in the PDB version 3.30 documentation at U{http://www.wwpdb.org/documentation/file-format/format33/sect4.html}. 391 392 393 @param lines: The lines of the PDB file excluding the sections prior to the hetrogen section. 394 @type lines: list of str 395 @return: The remaining PDB lines with the hetrogen records stripped. 396 @rtype: list of str 397 """ 398 399 # The ordered list of record names in the hetrogen section. 400 records = [ 401 'HET ', 402 'FORMUL', 403 'HETNAM', 404 'HETSYN' 405 ] 406 407 # Loop over the lines. 408 for i in range(len(lines)): 409 # No match, therefore assume to be out of the hetrogen section. 410 if lines[i][:6] not in records: 411 break 412 413 # Return the remaining lines. 414 return lines[i:]

415 416

417 - def _parse_pdb_misc(self, lines):

418 """Loop over and parse the PDB miscellaneous records. 419 420 These are the records identified in the PDB version 3.30 documentation at U{http://www.wwpdb.org/documentation/file-format/format33/sect7.html}. 421 422 423 @param lines: The lines of the PDB file excluding the sections prior to the miscellaneous section. 424 @type lines: list of str 425 @return: The remaining PDB lines with the miscellaneous records stripped. 426 @rtype: list of str 427 """ 428 429 # The ordered list of record names in the miscellaneous section. 430 records = [ 431 'SITE ' 432 ] 433 434 # Loop over the lines. 435 for i in range(len(lines)): 436 # No match, therefore assume to be out of the miscellaneous section. 437 if lines[i][:6] not in records: 438 break 439 440 # Return the remaining lines. 441 return lines[i:]

442 443

444 - def _parse_pdb_prim_struct(self, lines):

445 """Loop over and parse the PDB primary structure records. 446 447 These are the records identified in the PDB version 3.30 documentation at U{http://www.wwpdb.org/documentation/file-format/format33/sect3.html}. 448 449 450 @param lines: The lines of the PDB file excluding the title section. 451 @type lines: list of str 452 @return: The remaining PDB lines with the primary structure records stripped. 453 @rtype: list of str 454 """ 455 456 # The ordered list of record names in the primary structure section. 457 records = [ 458 'DBREF ', 459 'DBREF1', 460 'DBREF2', 461 'SEQADV', 462 'SEQRES', 463 'MODRES' 464 ] 465 466 # Loop over the lines. 467 for i in range(len(lines)): 468 # No match, therefore assume to be out of the primary structure section. 469 if lines[i][:6] not in records: 470 break 471 472 # Return the remaining lines. 473 return lines[i:]

474 475

476 - def _parse_pdb_ss(self, lines, read_mol=None, helices=None, sheets=None):

477 """Loop over and parse the PDB secondary structure records. 478 479 These are the records identified in the PDB version 3.30 documentation at U{http://www.wwpdb.org/documentation/file-format/format33/sect5.html}. 480 481 482 @param lines: The lines of the PDB file excluding the sections prior to the secondary structure section. 483 @type lines: list of str 484 @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. 485 @type read_mol: None, int, or list of int 486 @keyword helices: The list of helix data to append to. 487 @type helices: list 488 @keyword sheets: The list of sheet data to append to. 489 @type sheets: list 490 @return: The remaining PDB lines with the secondary structure records stripped. 491 @rtype: list of str 492 """ 493 494 # The ordered list of record names in the secondary structure section (the depreciated TURN record is also included to handle old PDB files). 495 records = [ 496 'HELIX ', 497 'SHEET ', 498 'TURN ' 499 ] 500 501 # The number of pre-existing molecules. 502 if not len(self.structural_data): 503 mol_num = 0 504 else: 505 mol_num = len(self.structural_data[0].mol) 506 507 # Loop over the lines. 508 for i in range(len(lines)): 509 # No match, therefore assume to be out of the secondary structure section. 510 if lines[i][:6] not in records: 511 break 512 513 # A helix. 514 if lines[i][:5] == 'HELIX': 515 # Parse the record. 516 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]) 517 518 # The molecule indices. 519 mol_init_index = self._pdb_chain_id_to_mol_index(init_chain_id) 520 mol_end_index = self._pdb_chain_id_to_mol_index(end_chain_id) 521 522 # Only load the desired molecule. 523 if read_mol != None: 524 if mol_init_index + 1 not in read_mol: 525 continue 526 if mol_end_index + 1 not in read_mol: 527 continue 528 529 # Store the data. 530 helices.append([helix_id, mol_init_index, init_res_name, init_seq_num, mol_end_index, end_res_name, end_seq_num, helix_class, length]) 531 532 # A sheet. 533 if lines[i][:5] == 'SHEET': 534 # Parse the record. 535 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]) 536 537 # The molecule indices. 538 mol_init_index = self._pdb_chain_id_to_mol_index(init_chain_id) 539 mol_end_index = self._pdb_chain_id_to_mol_index(end_chain_id) 540 541 # Only load the desired molecule. 542 if read_mol != None: 543 if mol_init_index + 1 not in read_mol: 544 continue 545 if mol_end_index + 1 not in read_mol: 546 continue 547 548 # New molecule indices based on currently loaded data. 549 mol_cur_index = None 550 if cur_chain_id: 551 mol_cur_index = self._pdb_chain_id_to_mol_index(cur_chain_id) 552 mol_prev_index = None 553 if prev_chain_id: 554 mol_prev_index = self._pdb_chain_id_to_mol_index(prev_chain_id) 555 556 # Store the data. 557 sheets.append([strand, sheet_id, num_strands, init_res_name, mol_init_index, init_seq_num, init_icode, end_res_name, mol_end_index, end_seq_num, end_icode, sense, cur_atom, cur_res_name, mol_cur_index, cur_res_seq, cur_icode, prev_atom, prev_res_name, mol_prev_index, prev_res_seq, prev_icode]) 558 559 # Return the remaining lines. 560 return lines[i:]

561 562

563 - def _parse_pdb_title(self, lines):

564 """Loop over and parse the PDB title records. 565 566 These are the records identified in the PDB version 3.30 documentation at U{http://www.wwpdb.org/documentation/file-format/format33/sect2.html}. 567 568 569 @param lines: All lines of the PDB file. 570 @type lines: list of str 571 @return: The remaining PDB lines with the title records stripped. 572 @rtype: list of str 573 """ 574 575 # The ordered list of (sometimes truncated) record names in the title section. 576 records = [ 577 'HEADER', 578 'OBSLTE', 579 'TITLE ', 580 'SPLT ', 581 'CAVEAT', 582 'COMPND', 583 'SOURCE', 584 'KEYWDS', 585 'EXPDTA', 586 'NUMMDL', 587 'MDLTYP', 588 'AUTHOR', 589 'REVDAT', 590 'SPRSDE', 591 'JRNL ', 592 'REMARK' 593 ] 594 595 # Loop over the lines. 596 for i in range(len(lines)): 597 # No match, therefore assume to be out of the title section. 598 if lines[i][:6] not in records: 599 break 600 601 # Return the remaining lines. 602 return lines[i:]

603 604

605 - def _parse_pdb_transform(self, lines):

606 """Loop over and parse the PDB transform records. 607 608 These are the records identified in the PDB version 3.30 documentation at U{http://www.wwpdb.org/documentation/file-format/format33/sect8.html}. 609 610 611 @param lines: The lines of the PDB file excluding the sections prior to the transform section. 612 @type lines: list of str 613 @return: The remaining PDB lines with the transform records stripped. 614 @rtype: list of str 615 """ 616 617 # The ordered list of record names in the transform section. 618 records = [ 619 'CRYST', 620 'MTRIX', 621 'ORIGX', 622 'SCALE', 623 ] 624 625 # Loop over the lines. 626 for i in range(len(lines)): 627 # No match, therefore assume to be out of the transform section. 628 if lines[i][0: 5] not in records: 629 break 630 631 # Return the remaining lines. 632 return lines[i:]

633 634

635 - def _parse_models_xyz(self, file_path, verbosity=1):

636 """Generator function for looping over the models in the XYZ file. 637 638 @param file_path: The full path of the XYZ file. 639 @type file_path: str 640 @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. 641 @type verbosity: int 642 @return: The model number and all the records for that model. 643 @rtype: tuple of int and array of str 644 """ 645 646 # Open the file. 647 file = open_read_file(file_path, verbosity=verbosity) 648 lines = file.readlines() 649 file.close() 650 651 # Check for empty files. 652 if lines == []: 653 raise RelaxError("The XYZ file is empty.") 654 655 # Init. 656 total_atom = 0 657 model = 0 658 records = [] 659 660 # Loop over the data. 661 for i in range(len(lines)): 662 num=0 663 word = lines[i].split() 664 # Find the total atom number and the first model. 665 if (i==0) and (len(word)==1): 666 try: 667 total_atom = int(word[0]) 668 num = 1 669 except: 670 raise RelaxError("The MODEL record " + repr(lines[i]) + " is corrupt, cannot read the XYZ file.") 671 672 # End of the model. 673 if (len(records) == total_atom): 674 # Yield the info 675 yield records 676 677 # Reset the records. 678 records = [] 679 680 # Skip all records prior to atom coordinates record. 681 if (len(word) != 4): 682 continue 683 684 # Append the line as a record of the model. 685 records.append(lines[i]) 686 687 # If records is not empty then there are no models, so yield the lot. 688 if len(records): 689 yield records

690 691

692 - def _parse_mols_pdb(self, records):

693 """Generator function for looping over the molecules in the PDB records of a model. 694 695 @param records: The list of PDB records for the model, or if no models exist the entire PDB file. 696 @type records: list of str 697 @return: The molecule number and all the records for that molecule. 698 @rtype: tuple of int and list of str 699 """ 700 701 # Check for empty records. 702 if records == []: 703 raise RelaxError("There are no PDB records for this model.") 704 705 # Init. 706 mol_count = 1 707 mol_records = [[]] 708 end = False 709 710 # Loop over the data. 711 for i in range(len(records)): 712 # A PDB termination record. 713 if records[i][:3] == 'END': 714 break 715 716 # A master record, so we are done. 717 if records[i][:6] == 'MASTER': 718 break 719 720 # A model termination record. 721 if records[i][:6] == 'ENDMDL': 722 end = True 723 724 # A molecule termination record with no trailing HETATM or CONECT. 725 elif i < len(records)-1 and records[i][:3] == 'TER' and not records[i+1][:6] == 'HETATM' and not records[i+1][:6] == 'CONECT': 726 end = True 727 728 # A HETATM followed by an ATOM record. 729 elif i < len(records)-1 and records[i][:6] == 'HETATM' and records[i+1][:4] == 'ATOM': 730 end = True 731 732 # End. 733 if end: 734 # Increment the molecule counter. 735 mol_count = mol_count + 1 736 737 # Reset the flag. 738 end = False 739 740 # The molecule number. 741 chain_id = records[i][21] 742 if chain_id == ' ': 743 mol_index = mol_count - 1 744 else: 745 mol_index = self._pdb_chain_id_to_mol_index(chain_id) 746 747 # Add a new records list as required. 748 while True: 749 if len(mol_records) <= mol_index: 750 mol_records.append([]) 751 else: 752 break 753 754 # Append the line as a record of the molecule. 755 mol_records[mol_index].append(records[i]) 756 757 # Loop over the molecules and yield the molecule number and records. 758 for i in range(len(mol_records)): 759 if mol_records[i] != []: 760 yield i+1, mol_records[i]

761 762

763 - def _pdb_chain_id_to_mol_index(self, chain_id=None):

764 """Convert the PDB chain ID into the molecule index in a regular way. 765 766 @keyword chain_id: The PDB chain ID string. 767 @type chain_id: str 768 @return: The corresponding molecule index. 769 @rtype: int 770 """ 771 772 # Initialise. 773 mol_index = 0 774 775 # Convert to the molecule index. 776 if chain_id: 777 mol_index = CHAIN_ID_LIST.index(chain_id) 778 779 # Return the index. 780 return mol_index

781 782

783 - def _residue_data(self, res_nums=None, res_names=None):

784 """Convert the residue info into a dictionary of unique residues with numbers as keys. 785 786 @keyword res_nums: The list of residue numbers. 787 @type res_nums: list of int 788 @keyword res_names: The list of residue names matching the numbers. 789 @type res_names: list of str 790 @return: The dictionary of residue names with residue numbers as keys. 791 @rtype: dict of str 792 """ 793 794 # Initialise. 795 data = {} 796 797 # Loop over the data. 798 for i in range(len(res_nums)): 799 # The residue data already exists. 800 if res_nums[i] in data: 801 continue 802 803 # Add the data. 804 data[res_nums[i]] = res_names[i] 805 806 # Return the dictionary. 807 return data

808 809

810 - def _validate_data_arrays(self, struct):

811 """Check the validity of the data arrays in the given structure object. 812 813 @param struct: The structural object. 814 @type struct: Structure_container instance 815 """ 816 817 # The number of atoms. 818 num = len(struct.atom_name) 819 820 # Check the other lengths. 821 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: 822 raise RelaxError("The structural data is invalid.")

823 824

825 - def _validate_records(self, lines):

826 """Make sure all PDB records are 80 char in length, padding with whitespace when needed. 827 828 All newline characters are stripped from the records as well. 829 830 831 @param lines: All lines of the PDB file. 832 @type lines: list of str 833 @return: The padded PDB lines. 834 @rtype: list of str 835 """ 836 837 # Loop over the lines. 838 for i in range(len(lines)): 839 # Strip the newline character. 840 lines[i] = lines[i].rstrip('\r\n') 841 842 # Pad if needed. 843 if len(lines[i]) != 80: 844 lines[i] = "%-80s" % lines[i] 845 846 # Return the fixed lines. 847 return lines

848 849

850 - def _mol_type(self, mol):

851 """Determine the type of molecule. 852 853 @param mol: The molecule data container. 854 @type mol: MolContainer instance 855 """ 856 857 # Amino acids. 858 aa = ['ALA', 'ARG', 'ASN', 'ASP', 'CYS', 'GLU', 'GLN', 'GLY', 'HIS', 'ILE', 'LEU', 'LYS', 'MET', 'PHE', 'PRO', 'SER', 'THR', 'TRP', 'TYR', 'VAL'] 859 860 # Set the molecule type to default to 'other'. 861 mol.type = 'other' 862 863 # Loop over the residues. 864 for res in mol.res_name: 865 # Protein. 866 if res in aa: 867 # Set the molecule type and return. 868 mol.type = 'protein' 869 return

870 871

872 - def _protein_connect(self, mol):

873 """Set up the connectivities for the protein. 874 875 @param mol: The molecule data container. 876 @type mol: MolContainer instance 877 """ 878 879 # Initialise some residue data. 880 curr_res_num = None 881 res_atoms = [] 882 883 # Loop over all atoms. 884 for i in range(len(mol.atom_num)): 885 # New residue. 886 if mol.res_num[i] != curr_res_num: 887 # Intra-residue connectivites. 888 if len(res_atoms): 889 self._protein_intra_connect(mol, res_atoms) 890 891 # Update the residue number. 892 curr_res_num = mol.res_num[i] 893 894 # Reset the residue atom index list. 895 res_atoms = [] 896 897 # Add the atom index to the list. 898 res_atoms.append(i) 899 900 # Last atom. 901 if i == len(mol.atom_num) - 1 and len(res_atoms): 902 self._protein_intra_connect(mol, res_atoms)

903 904

905 - def _protein_intra_connect(self, mol, res_atoms):

906 """Set up the connectivities for the protein. 907 908 @param mol: The molecule data container. 909 @type mol: MolContainer instance 910 @param res_atoms: The list of atom indices corresponding to the residue. 911 @type res_atoms: list of int 912 """ 913 914 # Back bond connectivity. 915 indices = { 916 'N': None, 917 'C': None, 918 'O': None, 919 'CA': None, 920 'HN': None, 921 'H': None, # Same as HN. 922 'HA': None 923 } 924 925 # Loop over all atoms to find the indices. 926 for index in res_atoms: 927 if mol.atom_name[index] in indices: 928 indices[mol.atom_name[index]] = index 929 930 # Connect the atom pairs. 931 pairs = [ 932 ['N', 'HN'], 933 ['N', 'H'], 934 ['N', 'CA'], 935 ['CA', 'HA'], 936 ['CA', 'C'], 937 ['C', 'O'] 938 ] 939 940 # Loop over the atoms pairs and connect them. 941 for pair in pairs: 942 if indices[pair[0]] != None and indices[pair[1]] != None: 943 mol.atom_connect(indices[pair[0]], indices[pair[1]])

944 945

946 - def _translate(self, data, format='str'):

947 """Convert the data into a format for writing to file. 948 949 @param data: The data to convert to the required format. 950 @type data: anything 951 @keyword format: The format to convert to. This can be 'str', 'float', or 'int'. 952 @type format: str 953 @return: The converted version of the data. 954 @rtype: str 955 """ 956 957 # Conversion to string. 958 if format == 'str': 959 # None values. 960 if data == None: 961 data = '' 962 963 # Force convert to string. 964 if not isinstance(data, str): 965 data = repr(data) 966 967 # Conversion to float. 968 if format == 'float': 969 # None values. 970 if data == None: 971 data = 0.0 972 973 # Force convert to float. 974 if not isinstance(data, float): 975 data = float(data) 976 977 # Return the converted data. 978 return data

979 980

981 - def _trim_helix(self, helix=None, trim_res_list=[], res_data=None):

982 """Trim the given helix based on the list of deleted residue numbers. 983 984 @keyword helix: The single helix metadata structure. 985 @type helix: list 986 @keyword trim_res_list: The list of residue numbers which no longer exist. 987 @type trim_res_list: list of int 988 @keyword res_data: The dictionary of residue names with residue numbers as keys. 989 @type res_data: dict of str 990 @return: The trimmed helix metadata structure, or None if the whole helix is to be deleted. 991 @rtype: list or None 992 """ 993 994 # Unpack the helix residue numbers. 995 start_res = helix[3] 996 end_res = helix[6] 997 998 # The reverse residue list. 999 trim_res_list_rev = deepcopy(trim_res_list) 1000 trim_res_list_rev.reverse() 1001 1002 # The helix residues. 1003 helix_res = list(range(start_res, end_res+1)) 1004 1005 # Trim forwards. 1006 for res_num in trim_res_list: 1007 if res_num == start_res: 1008 # Remove the residue. 1009 helix_res.pop(0) 1010 1011 # No helix left. 1012 if len(helix_res) == 0: 1013 break 1014 1015 # Realias the starting residue. 1016 start_res = helix_res[0] 1017 1018 # No helix left. 1019 if len(helix_res) == 0: 1020 return None 1021 1022 # Trim backwards. 1023 for res_num in trim_res_list_rev: 1024 if res_num == end_res: 1025 helix_res.pop(-1) 1026 end_res = helix_res[-1] 1027 1028 # Replace the starting and ending residues. 1029 if start_res != helix[3]: 1030 helix[3] = start_res 1031 helix[2] = res_data[start_res] 1032 if end_res != helix[6]: 1033 helix[6] = end_res 1034 helix[5] = res_data[end_res] 1035 1036 # The helix length. 1037 helix[-1] = len(helix_res) 1038 1039 # Return the modified helix. 1040 return helix

1041 1042

1043 - def _trim_sheet(self, sheet=None, trim_res_list=[], res_data=None):

1044 """Trim the given sheet based on the list of deleted residue numbers. 1045 1046 @keyword sheet: The single sheet metadata structure. 1047 @type sheet: list 1048 @keyword trim_res_list: The list of residue numbers which no longer exist. 1049 @type trim_res_list: list of int 1050 @keyword res_data: The dictionary of residue names with residue numbers as keys. 1051 @type res_data: dict of str 1052 @return: The trimmed sheet metadata structure, or None if the whole sheet is to be deleted. 1053 @rtype: list or None 1054 """ 1055 1056 # Unpack the sheet residue numbers. 1057 start_res = sheet[5] 1058 end_res = sheet[9] 1059 1060 # The reverse residue list. 1061 trim_res_list_rev = deepcopy(trim_res_list) 1062 trim_res_list_rev.reverse() 1063 1064 # The sheet residues. 1065 sheet_res = list(range(start_res, end_res+1)) 1066 1067 # Trim forwards. 1068 for res_num in trim_res_list: 1069 if res_num == start_res: 1070 # Remove the residue. 1071 sheet_res.pop(0) 1072 1073 # No sheet left. 1074 if len(sheet_res) == 0: 1075 break 1076 1077 # Realias the starting residue. 1078 start_res = sheet_res[0] 1079 1080 # No sheet left. 1081 if len(sheet_res) == 0: 1082 return None 1083 1084 # Trim backwards. 1085 for res_num in trim_res_list_rev: 1086 if res_num == end_res: 1087 sheet_res.pop(-1) 1088 end_res = sheet_res[-1] 1089 1090 # Replace the starting and ending residues. 1091 if start_res != sheet[5]: 1092 sheet[5] = start_res 1093 sheet[3] = res_data[start_res] 1094 if end_res != sheet[9]: 1095 sheet[9] = end_res 1096 sheet[7] = res_data[end_res] 1097 1098 # Return the modified sheet. 1099 return sheet

1100 1101

1102 - 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, sort=False):

1103 """Add a new atom to the structural data object. 1104 1105 @keyword mol_name: The name of the molecule. 1106 @type mol_name: str 1107 @keyword atom_name: The atom name, e.g. 'H1'. 1108 @type atom_name: str or None 1109 @keyword res_name: The residue name. 1110 @type res_name: str or None 1111 @keyword res_num: The residue number. 1112 @type res_num: int or None 1113 @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. 1114 @type pos: rank-1 or rank-2 array or list of float 1115 @keyword element: The element symbol. 1116 @type element: str or None 1117 @keyword atom_num: The atom number. 1118 @type atom_num: int or None 1119 @keyword chain_id: The chain identifier. 1120 @type chain_id: str or None 1121 @keyword segment_id: The segment identifier. 1122 @type segment_id: str or None 1123 @keyword pdb_record: The optional PDB record name, e.g. 'ATOM' or 'HETATM'. 1124 @type pdb_record: str or None 1125 @keyword sort: A flag which if True will cause the structural data to be sorted after adding the atom. 1126 @type sort: bool 1127 """ 1128 1129 # Add a model if not present. 1130 if len(self.structural_data) == 0: 1131 self.add_model() 1132 1133 # Check the position. 1134 if is_float(pos[0]): 1135 if len(pos) != 3: 1136 raise RelaxError("The single atomic position %s must be a 3D list." % pos) 1137 else: 1138 if len(pos) != len(self.structural_data): 1139 raise RelaxError("The %s atomic positions does not match the %s models present." % (len(pos), len(self.structural_data))) 1140 1141 # Loop over each model. 1142 for i in range(len(self.structural_data)): 1143 # Alias the model. 1144 model = self.structural_data[i] 1145 1146 # Specific molecule. 1147 mol = self.get_molecule(mol_name, model=model.num) 1148 1149 # Add the molecule, if it does not exist. 1150 if mol == None: 1151 self.add_molecule(name=mol_name) 1152 mol = self.get_molecule(mol_name, model=model.num) 1153 1154 # Split up the position if needed. 1155 if is_float(pos[0]): 1156 model_pos = pos 1157 else: 1158 model_pos = pos[i] 1159 1160 # Add the atom. 1161 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) 1162 1163 # Sort. 1164 if sort: 1165 mol._sort()

1166 1167

1168 - def add_coordinates(self, coord=None, mol_names=None, res_names=None, res_nums=None, atom_names=None, elements=None, set_mol_name=None, set_model_num=None):

1169 """Add a set of coordinates to the structural object. 1170 1171 @keyword coord: The array of atomic coordinates (first dimension is the atoms and the second are the coordinates). 1172 @type coord: numpy rank-2 float64 array 1173 @keyword mol_names: The list of molecule names corresponding to the first dimension of the coordinate array. 1174 @type mol_names: list of str 1175 @keyword res_names: The list of residue names corresponding to the first dimension of the coordinate array. 1176 @type res_names: list of str 1177 @keyword res_nums: The list of residue numbers corresponding to the first dimension of the coordinate array. 1178 @type res_nums: list of str 1179 @keyword atom_names: The list of atom names corresponding to the first dimension of the coordinate array. 1180 @type atom_names: list of str 1181 @keyword elements: The list of elements corresponding to the first dimension of the coordinate array. 1182 @type elements: list of str 1183 @keyword set_mol_name: Set the names of the molecules which are loaded. If set to None, then all molecule names must be identical and the new molecule will have the same name. 1184 @type set_mol_name: None or str 1185 @keyword set_model_num: Set the model number for the coordinates. 1186 @type set_model_num: None or int 1187 """ 1188 1189 # The new molecule name. 1190 if not set_mol_name: 1191 set_mol_name = mol_names[0] 1192 for name in mol_names: 1193 if set_mol_name != name: 1194 raise RelaxError("No unique molecule name can be found in the list %s for storing the new molecule structure.") 1195 1196 # No model set, so delete all current data. 1197 if set_mol_name == None and set_model_num == None: 1198 print("Deleting all structural data so that only the new molecule will be present.") 1199 self.delete(verbosity=0) 1200 1201 # Generate a molecule container for the new data. 1202 mol = MolContainer() 1203 1204 # Add the data. 1205 for i in range(len(coord)): 1206 mol.atom_add(atom_name=atom_names[i], res_name=res_names[i], res_num=res_nums[i], pos=coord[i], element=elements[i]) 1207 1208 # Create the structural data data structures. 1209 self.pack_structs([[mol]], orig_model_num=[None], set_model_num=[set_model_num], orig_mol_num=[[None]], set_mol_name=[set_mol_name])

1210 1211

1212 - def add_helix(self, res_start=None, res_end=None, mol_name=None):

1213 """Define new alpha helical secondary structure. 1214 1215 @keyword res_start: The residue number for the start of the helix. 1216 @type res_start: int 1217 @keyword res_end: The residue number for the end of the helix. 1218 @type res_end: int 1219 @keyword mol_name: Define the secondary structure for a specific molecule. 1220 @type mol_name: str or None 1221 """ 1222 1223 # Initialise the data structure if required. 1224 if not hasattr(self, 'helices'): 1225 self.helices = [] 1226 1227 # Use the first model. 1228 model = self.structural_data[0] 1229 1230 # Find the molecule to use. 1231 mol = None 1232 mol_name_list = [] 1233 for i in range(len(model.mol)): 1234 mol_name_list.append(model.mol[i].mol_name) 1235 if model.mol[i].mol_name == mol_name: 1236 mol = model.mol[i] 1237 mol_index = i 1238 if mol == None: 1239 raise RelaxError("Cannot find the molecule '%s' in %s." % (mol_name, mol_name_list)) 1240 1241 # Find the residue names. 1242 res_start_name = None 1243 res_end_name = None 1244 length = 0 1245 inside = False 1246 current_res_num = None 1247 for i in range(len(mol.res_num)): 1248 # The start of the helix. 1249 if mol.res_num[i] == res_start: 1250 res_start_name = mol.res_name[i] 1251 inside = True 1252 1253 # Residue count. 1254 if inside and current_res_num != mol.res_num[i]: 1255 current_res_num = mol.res_num[i] 1256 length += 1 1257 1258 # The end of the helix. 1259 if mol.res_num[i] == res_end: 1260 res_end_name = mol.res_name[i] 1261 inside = False 1262 1263 # Generate the list of [helix_id, mol_init_index, init_res_name, init_seq_num, mol_end_index, end_res_name, end_seq_num, helix_class, length]. 1264 helix = [len(self.helices)+1, mol_index, res_start_name, res_start, mol_index, res_end_name, res_end, 1, length] 1265 1266 # Add the data. 1267 self.helices.append(helix)

1268 1269

1270 - def add_model(self, model=None, coords_from=None):

1271 """Add a new model to the store. 1272 1273 The new model will be constructed 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. 1274 1275 @keyword model: The number of the model to create. 1276 @type model: int or None 1277 @keyword coords_from: The model number to take the coordinates from. 1278 @type coords_from: int or None 1279 @return: The model container. 1280 @rtype: ModelContainer instance 1281 """ 1282 1283 # Check if the model currently exists. 1284 if model != None: 1285 for i in range(len(self.structural_data)): 1286 if model == self.structural_data[i].num: 1287 raise RelaxError("The model '%s' already exists." % model) 1288 1289 # Set the model number for a single structure, if necessary. 1290 if len(self.structural_data) == 1 and self.structural_data[0].num == None: 1291 self.structural_data[0].num = coords_from 1292 self.structural_data.current_models[0] = coords_from 1293 1294 # Add a new model. 1295 model = self.structural_data.add_item(model_num=model) 1296 1297 # The model to duplicate. 1298 model_from = None 1299 if coords_from == None: 1300 model_from = self.structural_data[0] 1301 else: 1302 for i in range(len(self.structural_data)): 1303 if self.structural_data[i].num == coords_from: 1304 model_from = self.structural_data[i] 1305 break 1306 if not model_from: 1307 raise RelaxError("Cannot find model %i to copy the data from." % coords_from) 1308 1309 # Duplicate all data from the MolList object down. 1310 for mol_index in range(len(model_from.mol)): 1311 # Create a new molecule container. 1312 model.mol.add_item(mol_name=model_from.mol[mol_index].mol_name, mol_cont=MolContainer()) 1313 mol = model.mol[mol_index] 1314 mol_from = model_from.mol[mol_index] 1315 1316 # Loop over the atomic data. 1317 for i in range(len(mol_from.atom_num)): 1318 mol.atom_num.append(mol_from.atom_num[i]) 1319 mol.atom_name.append(mol_from.atom_name[i]) 1320 mol.bonded.append(mol_from.bonded[i]) 1321 mol.chain_id.append(mol_from.chain_id[i]) 1322 mol.element.append(mol_from.element[i]) 1323 mol.pdb_record.append(mol_from.pdb_record[i]) 1324 mol.res_name.append(mol_from.res_name[i]) 1325 mol.res_num.append(mol_from.res_num[i]) 1326 mol.seg_id.append(mol_from.seg_id[i]) 1327 mol.x.append(mol_from.x[i]) 1328 mol.y.append(mol_from.y[i]) 1329 mol.z.append(mol_from.z[i]) 1330 1331 # Return the model. 1332 return self.structural_data[-1]

1333 1334

1335 - def add_molecule(self, model_num=None, name=None):

1336 """Add a new molecule to the store. 1337 1338 @keyword model_num: The optional model to add the molecule to. If not supplied, the molecule will be added to all models. 1339 @type model_num: None or int 1340 @keyword name: The molecule identifier string. 1341 @type name: str 1342 """ 1343 1344 # Add a model if necessary. 1345 if len(self.structural_data) == 0: 1346 self.add_model() 1347 1348 # Add the molecule to each model. 1349 for model in self.model_loop(model=model_num): 1350 model.mol.add_item(mol_name=name, mol_cont=MolContainer())

1351 1352

1353 - def add_sheet(self, strand=1, sheet_id='A', strand_count=2, strand_sense=0, res_start=None, res_end=None, mol_name=None, current_atom=None, prev_atom=None):

1354 """Define new alpha helical secondary structure. 1355 1356 @keyword strand: Strand number which starts at 1 for each strand within a sheet and increases by one. 1357 @type strand: int 1358 @keyword sheet_id: The sheet identifier. To match the PDB standard, sheet IDs should range from 'A' to 'Z'. 1359 @type sheet_id: str 1360 @keyword strand_count: The number of strands in the sheet. 1361 @type strand_count: int 1362 @keyword strand_sense: Sense of strand with respect to previous strand in the sheet. 0 if first strand, 1 if parallel, -1 if anti-parallel. 1363 @type strand_sense: int 1364 @keyword res_start: The residue number for the start of the sheet. 1365 @type res_start: int 1366 @keyword res_end: The residue number for the end of the sheet. 1367 @type res_end: int 1368 @keyword mol_name: Define the secondary structure for a specific molecule. 1369 @type mol_name: str or None 1370 @keyword current_atom: The name of the first atom in the current strand, to link the current back to the previous strand. 1371 @type current_atom: str or None 1372 @keyword prev_atom: The name of the last atom in the previous strand, to link the current back to the previous strand. 1373 @type prev_atom: str or None 1374 """ 1375 1376 # Initialise the data structure if required. 1377 if not hasattr(self, 'sheets'): 1378 self.sheets = [] 1379 1380 # Use the first model. 1381 model = self.structural_data[0] 1382 1383 # Find the molecule to use. 1384 mol = None 1385 mol_name_list = [] 1386 for i in range(len(model.mol)): 1387 mol_name_list.append(model.mol[i].mol_name) 1388 if model.mol[i].mol_name == mol_name: 1389 mol = model.mol[i] 1390 mol_index = i 1391 if mol == None: 1392 raise RelaxError("Cannot find the molecule '%s' in %s." % (mol_name, mol_name_list)) 1393 1394 # Find the residue names. 1395 res_start_name = None 1396 res_end_name = None 1397 length = 0 1398 inside = False 1399 current_res_num = None 1400 for i in range(len(mol.res_num)): 1401 # The start of the sheet. 1402 if mol.res_num[i] == res_start: 1403 res_start_name = mol.res_name[i] 1404 inside = True 1405 1406 # Residue count. 1407 if inside and current_res_num != mol.res_num[i]: 1408 current_res_num = mol.res_num[i] 1409 length += 1 1410 1411 # The end of the sheet. 1412 if mol.res_num[i] == res_end: 1413 res_end_name = mol.res_name[i] 1414 inside = False 1415 1416 # Current strand info. 1417 num_strands = strand_count 1418 init_res_name = res_start_name 1419 mol_init_index = mol_index 1420 init_seq_num = res_start 1421 init_icode = None 1422 end_res_name = res_end_name 1423 mol_end_index = mol_index 1424 end_seq_num = res_end 1425 end_icode = None 1426 sense = strand_sense 1427 1428 # Current sheet info. 1429 cur_atom_name = None 1430 cur_res_name = None 1431 mol_cur_index = None 1432 cur_res_seq = None 1433 cur_icode = None 1434 prev_atom_name = None 1435 prev_res_name = None 1436 mol_prev_index = None 1437 prev_res_seq = None 1438 prev_icode = None 1439 if strand > 1: 1440 # The current strand. 1441 cur_atom_name = current_atom 1442 cur_res_name = res_start_name 1443 mol_cur_index = mol_index 1444 cur_res_seq = res_start 1445 cur_icode = None 1446 1447 # The previous strand. 1448 for i in reversed(range(len(self.sheets))): 1449 if self.sheets[i][1] == sheet_id: 1450 prev_atom_name = prev_atom 1451 prev_res_name = self.sheets[i][7] 1452 mol_prev_index = self.sheets[i][4] 1453 prev_res_seq = self.sheets[i][9] 1454 prev_icode = None 1455 1456 1457 # Generate the list of [strand, sheet_id, num_strands, init_res_name, mol_init_index, init_seq_num, init_icode, end_res_name, mol_end_index, end_seq_num, end_icode, sense, cur_atom, cur_res_name, mol_cur_index, cur_res_seq, cur_icode, prev_atom_name, prev_res_name, mol_prev_index, prev_res_seq, prev_icode]. 1458 sheet = [strand, sheet_id, num_strands, init_res_name, mol_init_index, init_seq_num, init_icode, end_res_name, mol_end_index, end_seq_num, end_icode, sense, cur_atom_name, cur_res_name, mol_cur_index, cur_res_seq, cur_icode, prev_atom_name, prev_res_name, mol_prev_index, prev_res_seq, prev_icode] 1459 1460 # Add the data. 1461 self.sheets.append(sheet)

1462 1463

1464 - def are_bonded(self, atom_id1=None, atom_id2=None):

1465 """Determine if two atoms are directly bonded to each other. 1466 1467 @keyword atom_id1: The molecule, residue, and atom identifier string of the first atom. 1468 @type atom_id1: str 1469 @keyword atom_id2: The molecule, residue, and atom identifier string of the second atom. 1470 @type atom_id2: str 1471 @return: True if the atoms are directly bonded. 1472 @rtype: bool 1473 """ 1474 1475 # Generate the selection objects. 1476 sel_obj1 = Selection(atom_id1) 1477 sel_obj2 = Selection(atom_id2) 1478 1479 # Build the connectivities if needed. 1480 for mol in self.structural_data[0].mol: 1481 for i in range(len(mol.atom_num)): 1482 if not len(mol.bonded[i]): 1483 self._find_bonded_atoms(i, mol, radius=2) 1484 1485 # Loop over the molecules. 1486 for mol in self.structural_data[0].mol: 1487 # Skip non-matching molecules. 1488 if not sel_obj1.contains_mol(mol.mol_name): 1489 continue 1490 if not sel_obj2.contains_mol(mol.mol_name): 1491 continue 1492 1493 # Find the first atom. 1494 index1 = None 1495 index2 = None 1496 for i in range(len(mol.atom_num)): 1497 # Matching first atom. 1498 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): 1499 index1 = i 1500 1501 # Matching second atom. 1502 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): 1503 index2 = i 1504 1505 # Nothing left to do. 1506 if index1 != None and index2 != None: 1507 break 1508 1509 # Connectivities exist. 1510 if index1 < len(mol.bonded): 1511 if index2 in mol.bonded[index1]: 1512 return True 1513 else: 1514 return False

1515 1516

1517 - def are_bonded_index(self, mol_index1=None, atom_index1=None, mol_index2=None, atom_index2=None):

1518 """Determine if two atoms, given as indices, are directly bonded to each other. 1519 1520 @keyword mol_index1: The molecule index of the first atom. 1521 @type mol_index1: int 1522 @keyword atom_index1: The index of the first atom. 1523 @type atom_index1: int 1524 @keyword mol_index2: The molecule index of the second atom. 1525 @type mol_index2: int 1526 @keyword atom_index2: The index of the second atom. 1527 @type atom_index2: int 1528 @return: True if the atoms are directly bonded. 1529 @rtype: bool 1530 """ 1531 1532 # Alias the molecule. 1533 mol1 = self.structural_data[0].mol[mol_index1] 1534 mol2 = self.structural_data[0].mol[mol_index2] 1535 1536 # Build the connectivities if needed. 1537 if not len(mol1.bonded[atom_index1]): 1538 self._find_bonded_atoms(atom_index1, mol1, radius=2) 1539 if not len(mol2.bonded[atom_index2]): 1540 self._find_bonded_atoms(atom_index2, mol2, radius=2) 1541 1542 # Is the second atom in the bonded list of the first? 1543 if atom_index2 in mol1.bonded[atom_index1]: 1544 return True 1545 1546 # No! 1547 return False

1548 1549

1550 - def atom_loop(self, selection=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):

1551 """Generator function for looping over all atoms in the internal relax structural object. 1552 1553 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: 1554 1555 1. Model number. 1556 2. Molecule name. 1557 3. Residue number. 1558 4. Residue name. 1559 5. Atom number. 1560 6. Atom name. 1561 7. The element name (its atomic symbol and optionally the isotope, e.g. 'N', 'Mg', '17O', '13C', etc). 1562 8. The position of the atom in Euclidean space. 1563 1564 1565 @keyword selection: The internal structural selection object. This is obtained by calling the selection() method with the atom ID string. 1566 @type selection: lib.structure.internal.Internal_selection instance 1567 @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. 1568 @type str_id: str, int, or None 1569 @keyword model_num: Only loop over a specific model. 1570 @type model_num: int or None 1571 @keyword mol_name_flag: A flag which if True will cause the molecule name to be yielded. 1572 @type mol_name_flag: bool 1573 @keyword res_num_flag: A flag which if True will cause the residue number to be yielded. 1574 @type res_num_flag: bool 1575 @keyword res_name_flag: A flag which if True will cause the residue name to be yielded. 1576 @type res_name_flag: bool 1577 @keyword atom_num_flag: A flag which if True will cause the atom number to be yielded. 1578 @type atom_num_flag: bool 1579 @keyword atom_name_flag: A flag which if True will cause the atom name to be yielded. 1580 @type atom_name_flag: bool 1581 @keyword element_flag: A flag which if True will cause the element name to be yielded. 1582 @type element_flag: bool 1583 @keyword pos_flag: A flag which if True will cause the atomic position to be yielded. 1584 @type pos_flag: bool 1585 @keyword mol_index_flag: A flag which if True will cause the molecule index to be yielded. 1586 @type mol_index_flag: bool 1587 @keyword index_flag: A flag which if True will cause the atomic index to be yielded. 1588 @type index_flag: bool 1589 @keyword ave: A flag which if True will result in this method returning the average atom properties across all loaded structures. 1590 @type ave: bool 1591 @return: A tuple of atomic information, as described in the docstring. 1592 @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). 1593 """ 1594 1595 # Check that the structure is loaded. 1596 if not len(self.structural_data): 1597 raise RelaxNoPdbError 1598 1599 # Obtain all data from the first model (except the position data). 1600 model = self.structural_data[0] 1601 1602 # Loop over all molecules and atoms in the selection. 1603 for mol_index, i in selection.loop(): 1604 mol = model.mol[mol_index] 1605 1606 # Initialise. 1607 res_num = mol.res_num[i] 1608 res_name = mol.res_name[i] 1609 atom_num = mol.atom_num[i] 1610 atom_name = mol.atom_name[i] 1611 element = mol.element[i] 1612 1613 # The atom position. 1614 if pos_flag: 1615 # Average the position. 1616 if ave: 1617 # Initialise. 1618 pos = zeros(3, float64) 1619 1620 # Loop over the models. 1621 for model in self.model_loop(model=model_num): 1622 # Alias. 1623 mol2 = model.mol[mol_index] 1624 1625 # Some sanity checks. 1626 if mol2.atom_num[i] != atom_num: 1627 raise RelaxError("The loaded structures do not contain the same atoms. The average structural properties can not be calculated.") 1628 1629 # Sum the atom positions. 1630 pos = pos + array([mol2.x[i], mol2.y[i], mol2.z[i]], float64) 1631 1632 # Average the position array (divide by the number of models). 1633 pos = pos / len(self.structural_data) 1634 1635 # All positions. 1636 else: 1637 # Initialise. 1638 pos = [] 1639 1640 # Loop over the models. 1641 for model in self.model_loop(model=model_num): 1642 # Alias. 1643 mol2 = model.mol[mol_index] 1644 1645 # Append the position. 1646 pos.append([mol2.x[i], mol2.y[i], mol2.z[i]]) 1647 1648 # Convert. 1649 pos = array(pos, float64) 1650 1651 # The molecule name. 1652 mol_name = mol.mol_name 1653 1654 # Build the tuple to be yielded. 1655 atomic_tuple = () 1656 if mol_name_flag: 1657 atomic_tuple = atomic_tuple + (mol_name,) 1658 if res_num_flag: 1659 atomic_tuple = atomic_tuple + (res_num,) 1660 if res_name_flag: 1661 atomic_tuple = atomic_tuple + (res_name,) 1662 if atom_num_flag: 1663 atomic_tuple = atomic_tuple + (atom_num,) 1664 if atom_name_flag: 1665 atomic_tuple = atomic_tuple + (atom_name,) 1666 if element_flag: 1667 atomic_tuple = atomic_tuple + (element,) 1668 if pos_flag: 1669 atomic_tuple = atomic_tuple + (pos,) 1670 if mol_index_flag: 1671 atomic_tuple += (mol_index,) 1672 if index_flag: 1673 atomic_tuple += (i,) 1674 1675 # Yield the information. 1676 if len(atomic_tuple) == 1: 1677 atomic_tuple = atomic_tuple[0] 1678 yield atomic_tuple

1679 1680

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

1682 """Find the bond vector between the atoms of 'attached_atom' and 'atom_id'. 1683 1684 @keyword attached_atom: The name of the bonded atom. 1685 @type attached_atom: str 1686 @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. 1687 @type model_num: None or int 1688 @keyword mol_name: The name of the molecule that attached_atom belongs to. 1689 @type mol_name: str 1690 @keyword res_num: The number of the residue that attached_atom belongs to. 1691 @type res_num: str 1692 @keyword res_name: The name of the residue that attached_atom belongs to. 1693 @type res_name: str 1694 @keyword spin_num: The number of the spin that attached_atom is attached to. 1695 @type spin_num: str 1696 @keyword spin_name: The name of the spin that attached_atom is attached to. 1697 @type spin_name: str 1698 @keyword return_name: A flag which if True will cause the name of the attached atom to be returned together with the bond vectors. 1699 @type return_name: bool 1700 @keyword return_warnings: A flag which if True will cause warning messages to be returned. 1701 @type return_warnings: bool 1702 @return: The list of bond vectors for each model. 1703 @rtype: list of numpy arrays (or a tuple if return_name or return_warnings are set) 1704 """ 1705 1706 # Initialise some objects. 1707 vectors = [] 1708 attached_name = None 1709 warnings = None 1710 1711 # Use the first model for the atom matching. 1712 model = self.structural_data[0] 1713 1714 # Loop over the molecules. 1715 for mol_index in range(len(model.mol)): 1716 # Alias. 1717 mol = model.mol[mol_index] 1718 1719 # Skip non-matching molecules. 1720 if mol_name and mol_name != mol.mol_name: 1721 continue 1722 1723 # Find the atomic index of the base atom. 1724 index = None 1725 for i in range(len(mol.atom_name)): 1726 # Residues don't match. 1727 if (res_num != None and mol.res_num[i] != res_num) or (res_name != None and mol.res_name[i] != res_name): 1728 continue 1729 1730 # Atoms don't match. 1731 if (spin_num != None and mol.atom_num[i] != spin_num) or (spin_name != None and mol.atom_name[i] != spin_name): 1732 continue 1733 1734 # Update the index and stop searching. 1735 index = i 1736 break 1737 1738 # Found the atom. 1739 if index != None: 1740 # Loop over the models. 1741 for model in self.model_loop(model=model_num): 1742 # Alias. 1743 mol = model.mol[mol_index] 1744 1745 # Get the atom bonded to this model/molecule/residue/atom. 1746 bonded_num, bonded_name, element, pos, attached_name, warnings = self._bonded_atom(attached_atom, index, mol) 1747 1748 # No bonded atom. 1749 if (bonded_num, bonded_name, element) == (None, None, None): 1750 continue 1751 1752 # The bond vector. 1753 vector = array(pos, float64) - array([mol.x[index], mol.y[index], mol.z[index]], float64) 1754 1755 # Append the vector to the vectors array. 1756 vectors.append(vector) 1757 1758 # Not found. 1759 else: 1760 warnings = "Cannot find the atom in the structure" 1761 1762 # Build the tuple to be yielded. 1763 data = (vectors,) 1764 if return_name: 1765 data = data + (attached_name,) 1766 if return_warnings: 1767 data = data + (warnings,) 1768 1769 # Return the data. 1770 return data

1771 1772

1773 - def collapse_ensemble(self, model_num=None, model_to=1):

1774 """Collapse the ensemble into a single model. 1775 1776 @keyword model_num: The number of the model to keep. All other models will be removed. 1777 @type model_num: int 1778 @keyword model_to: The model number for the sole remaining model. 1779 @type model_to: int 1780 """ 1781 1782 # Store all the model numbers. 1783 models = [] 1784 for model_cont in self.model_loop(): 1785 if model_cont.num != model_num: 1786 models.append(model_cont.num) 1787 1788 # Delete all models. 1789 for model in models: 1790 self.delete(model) 1791 1792 # Renumber the remaining model. 1793 self.set_model(model_orig=model_num, model_new=model_to)

1794 1795

1796 - def connect_atom(self, mol_name=None, index1=None, index2=None):

1797 """Connect two atoms in the structural data object. 1798 1799 @keyword mol_name: The name of the molecule. 1800 @type mol_name: str 1801 @keyword index1: The global index of the first atom. 1802 @type index1: str 1803 @keyword index2: The global index of the first atom. 1804 @type index2: str 1805 """ 1806 1807 # Add the molecule, if it does not exist. 1808 if self.get_molecule(mol_name) == None: 1809 self.add_molecule(name=mol_name) 1810 1811 # Loop over each model. 1812 for model in self.structural_data: 1813 # Specific molecule. 1814 mol = self.get_molecule(mol_name) 1815 1816 # Add the atom. 1817 mol.atom_connect(index1=index1, index2=index2)

1818 1819

1820 - def delete(self, model=None, selection=None, verbosity=1):

1821 """Deletion of structural information. 1822 1823 @keyword model: Individual structural models from a loaded ensemble can be deleted by specifying the model number. 1824 @type model: None or int 1825 @keyword selection: The internal structural selection object. This is obtained by calling the selection() method with the atom ID string. 1826 @type selection: lib.structure.internal.Internal_selection instance 1827 @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. 1828 @type verbosity: int 1829 """ 1830 1831 # All data. 1832 if model == None and selection == None: 1833 # Printout. 1834 if verbosity: 1835 print("Deleting the following structural data:\n") 1836 print(self.structural_data) 1837 1838 # Delete the structural data. 1839 del self.structural_data 1840 1841 # Initialise the empty model list. 1842 self.structural_data = ModelList() 1843 1844 # Delete a whole model. 1845 elif selection == None: 1846 self.structural_data.delete_model(model_num=model) 1847 1848 # Atom subset deletion. 1849 else: 1850 # Loop over the atoms and find the indices of the atoms to delete. 1851 indices = {} 1852 for mol_index, i in selection.loop(): 1853 if mol_index not in indices: 1854 indices[mol_index] = [] 1855 indices[mol_index].append(i) 1856 for mol_index in indices: 1857 indices[mol_index].reverse() 1858 1859 # Loop over the models. 1860 del_res_nums = [] 1861 for model_cont in self.model_loop(): 1862 # Skip models. 1863 if model != None and model_cont.num == model: 1864 continue 1865 1866 # Loop over the molecules. 1867 for mol_index in indices: 1868 mol = model_cont.mol[mol_index] 1869 1870 # Generate a residue data dictionary for the metadata trimming (prior to atom deletion). 1871 res_data = self._residue_data(res_nums=mol.res_num, res_names=mol.res_name) 1872 1873 # Loop over the reverse indices and pop out the data. 1874 for i in indices[mol_index]: 1875 mol.atom_num.pop(i) 1876 mol.atom_name.pop(i) 1877 mol.bonded.pop(i) 1878 mol.chain_id.pop(i) 1879 mol.element.pop(i) 1880 mol.pdb_record.pop(i) 1881 mol.res_name.pop(i) 1882 res_num = mol.res_num.pop(i) 1883 mol.seg_id.pop(i) 1884 mol.x.pop(i) 1885 mol.y.pop(i) 1886 mol.z.pop(i) 1887 1888 # The residue no longer exists. 1889 if res_num not in mol.res_num and res_num not in del_res_nums: 1890 del_res_nums.append(res_num) 1891 1892 # Second atom of the bonded pair. 1893 for j in range(len(mol.bonded)): 1894 if i in mol.bonded[j]: 1895 mol.bonded[j].pop(mol.bonded[j].index(i)) 1896 1897 # Update the bonded lists, as the indices need to be shifted. 1898 for j in range(i, len(mol.bonded)): 1899 for k in range(len(mol.bonded[j])): 1900 mol.bonded[j][k] -= 1 1901 1902 # Reset the metadata if nothing remains. 1903 if mol.atom_num == []: 1904 if hasattr(mol, 'file_name'): 1905 del mol.file_name 1906 if hasattr(mol, 'file_path'): 1907 del mol.file_path 1908 if hasattr(mol, 'file_mol_num'): 1909 del mol.file_mol_num 1910 if hasattr(mol, 'file_model'): 1911 del mol.file_model 1912 1913 # Nothing more to do. 1914 if not len(del_res_nums): 1915 return 1916 if model != None and len(self.structural_data) > 1: 1917 return 1918 1919 # Fix the deleted residue number order. 1920 del_res_nums.reverse() 1921 1922 # Handle the helix metadata. 1923 if hasattr(self, 'helices'): 1924 del_helix_indices = [] 1925 for i in range(len(self.helices)): 1926 # Trim the helix. 1927 helix = self._trim_helix(helix=self.helices[i], trim_res_list=del_res_nums, res_data=res_data) 1928 1929 # Trimmed helix. 1930 if helix != None: 1931 self.helices[i] = helix 1932 1933 # No helix left. 1934 else: 1935 del_helix_indices.append(i) 1936 1937 # Loop over the reverse helix indices and pop out the data. 1938 del_helix_indices.reverse() 1939 for i in del_helix_indices: 1940 self.helices.pop(i) 1941 1942 # Handle the sheet metadata. 1943 if hasattr(self, 'sheets'): 1944 del_sheet_indices = [] 1945 for i in range(len(self.sheets)): 1946 # Trim the sheet. 1947 sheet = self._trim_sheet(sheet=self.sheets[i], trim_res_list=del_res_nums, res_data=res_data) 1948 1949 # Trimmed sheet. 1950 if sheet != None: 1951 self.sheets[i] = sheet 1952 1953 # No sheet left. 1954 else: 1955 del_sheet_indices.append(i) 1956 1957 # Loop over the reverse sheet indices and pop out the data. 1958 del_sheet_indices.reverse() 1959 for i in del_sheet_indices: 1960 self.sheets.pop(i)

1961 1962

1963 - def delete_ss(self, verbosity=1):

1964 """Deletion of all secondary structure information. 1965 1966 @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. 1967 @type verbosity: int 1968 """ 1969 1970 # Reset the structures. 1971 self.helices = [] 1972 self.sheets = []

1973 1974

1975 - def empty(self):

1976 """Report if the structural data structure is empty or not. 1977 1978 @return: True if empty, False otherwise. 1979 @rtype: bool 1980 """ 1981 1982 # Check the ModelList structure. 1983 if len(self.structural_data) == 0: 1984 return True 1985 else: 1986 return False

1987 1988

1989 - def from_xml(self, str_node, dir=None, file_version=1):

1990 """Recreate the structural object from the XML structural object node. 1991 1992 @param str_node: The structural object XML node. 1993 @type str_node: xml.dom.minicompat.Element instance 1994 @keyword dir: The name of the directory containing the results file. 1995 @type dir: str 1996 @keyword file_version: The relax XML version of the XML file. 1997 @type file_version: int 1998 """ 1999 2000 # Recreate all base objects (i.e. metadata). 2001 xml_to_object(str_node, self, file_version=file_version, blacklist=['model', 'displacements']) 2002 2003 # Recreate the model / molecule data structure. 2004 model_nodes = str_node.getElementsByTagName('model') 2005 self.structural_data.from_xml(model_nodes, file_version=file_version) 2006 2007 # The displacement structure. 2008 disp_nodes = str_node.getElementsByTagName('displacements') 2009 if len(disp_nodes): 2010 # Initialise the object. 2011 self.displacements = Displacements() 2012 2013 # Recreate the molecule data structures for the current model. 2014 self.displacements.from_xml(disp_nodes[0], file_version=file_version)

2015 2016

2017 - def get_model(self, model):

2018 """Return or create the model. 2019 2020 @param model: The model number. 2021 @type model: int or None 2022 @return: The ModelContainer corresponding to the model number or that newly created. 2023 @rtype: ModelContainer instance 2024 """ 2025 2026 # Check if the target is a single model. 2027 if model == None and self.num_models() > 1: 2028 raise RelaxError("The target model cannot be determined as there are %s models already present." % self.num_modes()) 2029 2030 # No model specified. 2031 if model == None: 2032 # Create the first model, if necessary. 2033 if self.num_models(): 2034 self.structural_data.add_item() 2035 2036 # Alias the first model. 2037 model_cont = self.structural_data[0] 2038 2039 # The model has been specified. 2040 else: 2041 # Get the preexisting model. 2042 found = False 2043 for model_cont in self.structural_data: 2044 if model_cont.num == model: 2045 found = True 2046 break 2047 2048 # Add the model if it doesn't exist. 2049 if not found: 2050 self.structural_data.add_item(model) 2051 model_cont = self.structural_data[-1] 2052 2053 # Return the container. 2054 return model_cont

2055 2056

2057 - def get_molecule(self, molecule, model=None):

2058 """Return the molecule. 2059 2060 Only one model can be specified. 2061 2062 2063 @param molecule: The molecule name. 2064 @type molecule: int or None 2065 @keyword model: The model number. 2066 @type model: int or None 2067 @raises RelaxError: If the model is not specified and there is more than one model loaded. 2068 @return: The MolContainer corresponding to the molecule name and model number. 2069 @rtype: MolContainer instance or None 2070 """ 2071 2072 # Check if the target is a single molecule. 2073 if model == None and self.num_models() > 1: 2074 raise RelaxError("The target molecule cannot be determined as there are %s models already present." % self.num_models()) 2075 2076 # Check the model argument. 2077 if not isinstance(model, int) and not model == None: 2078 raise RelaxNoneIntError 2079 2080 # No models. 2081 if not len(self.structural_data): 2082 return 2083 2084 # Loop over the models. 2085 for model_cont in self.model_loop(model): 2086 # Loop over the molecules. 2087 for mol in model_cont.mol: 2088 # Return the matching molecule. 2089 if mol.mol_name == molecule: 2090 return mol

2091 2092

2093 - def has_molecule(self, model_num=None, name=None):

2094 """Check if the molecule name exists. 2095 2096 @keyword model_num: The optional model to check. If not supplied, the molecule will be searched for across all models. 2097 @type model_num: None or int 2098 @param name: The molecule name. 2099 @type name: str 2100 @return: True if the molecule exists, False otherwise. 2101 @rtype: bool 2102 """ 2103 2104 # No models. 2105 if not len(self.structural_data): 2106 return False 2107 2108 # Loop over the models. 2109 for model_cont in self.model_loop(model=model_num): 2110 # Loop over the molecules. 2111 for mol in model_cont.mol: 2112 # Matching molecule. 2113 if mol.mol_name == name: 2114 return True 2115 2116 # No match. 2117 return False

2118 2119

2120 - def load_gaussian(self, file_path, set_mol_name=None, set_model_num=None, verbosity=False):

2121 """Method for loading structures from a Gaussian log file. 2122 2123 @param file_path: The full path of the Gaussian log file. 2124 @type file_path: str 2125 @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. 2126 @type set_mol_name: None, str, or list of str 2127 @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. 2128 @type set_model_num: None, int, or list of int 2129 @keyword verbosity: A flag which if True will cause messages to be printed. 2130 @type verbosity: bool 2131 @return: The status of the loading of the Gaussian log file. 2132 @rtype: bool 2133 """ 2134 2135 # Initial printout. 2136 if verbosity: 2137 print("\nInternal relax Gaussian log parser.") 2138 2139 # Test if the file exists. 2140 if not access(file_path, F_OK): 2141 # Exit indicating failure. 2142 return False 2143 2144 # Separate the file name and path. 2145 path, file_name = os.path.split(file_path) 2146 2147 # The absolute path. 2148 path_abs = abspath(curdir) + sep + path 2149 2150 # Set up the molecule name data structure. 2151 if set_mol_name: 2152 if not isinstance(set_mol_name, list): 2153 set_mol_name = [set_mol_name] 2154 else: 2155 set_mol_name = [file_root(file_name) + '_mol1'] 2156 2157 # Set up the model number data structure. 2158 if set_model_num: 2159 if not isinstance(set_model_num, list): 2160 set_model_num = [set_model_num] 2161 else: 2162 set_model_num = [1] 2163 2164 # Loop over all models in the Gaussian log file, doing nothing so the last model records are stored. 2165 for model_records in self._parse_models_gaussian(file_path, verbosity=verbosity): 2166 pass 2167 2168 # Generate the molecule container. 2169 mol = MolContainer() 2170 2171 # Fill the molecular data object. 2172 mol.fill_object_from_gaussian(model_records) 2173 2174 # Create the structural data data structures. 2175 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) 2176 2177 # Loading worked. 2178 return True

2179 2180

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

2182 """Method for loading structures from a PDB file. 2183 2184 @param file_path: The full path of the PDB file. 2185 @type file_path: str 2186 @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. 2187 @type read_mol: None, int, or list of int 2188 @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. 2189 @type set_mol_name: None, str, or list of str 2190 @keyword read_model: The PDB model to extract from the file. If set to None, then all models will be loaded. 2191 @type read_model: None, int, or list of int 2192 @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. 2193 @type set_model_num: None, int, or list of int 2194 @keyword alt_loc: The PDB ATOM record 'Alternate location indicator' field value to select which coordinates to use. 2195 @type alt_loc: str or None 2196 @keyword verbosity: A flag which if True will cause messages to be printed. 2197 @type verbosity: bool 2198 @keyword merge: A flag which if set to True will try to merge the PDB structure into the currently loaded structures. 2199 @type merge: bool 2200 @return: The status of the loading of the PDB file. 2201 @rtype: bool 2202 """ 2203 2204 # Initial printout. 2205 if verbosity: 2206 print("\nInternal relax PDB parser.") 2207 2208 # Test if the file exists. 2209 if not access(file_path, F_OK): 2210 # Exit indicating failure. 2211 return False 2212 2213 # Separate the file name and path. 2214 path, file = os.path.split(file_path) 2215 2216 # The absolute path. 2217 path_abs = abspath(curdir) + sep + path 2218 2219 # Convert the structure reading args into lists. 2220 if read_mol and not isinstance(read_mol, list): 2221 read_mol = [read_mol] 2222 if set_mol_name and not isinstance(set_mol_name, list): 2223 set_mol_name = [set_mol_name] 2224 if read_mol: 2225 set_mol_name *= len(read_mol) 2226 if read_model and not isinstance(read_model, list): 2227 read_model = [read_model] 2228 if set_model_num and not isinstance(set_model_num, list): 2229 set_model_num = [set_model_num] 2230 if read_model: 2231 set_model_num *= len(read_model) 2232 2233 # Open the PDB file. 2234 pdb_file = open_read_file(file_path, verbosity=verbosity) 2235 pdb_lines = pdb_file.readlines() 2236 pdb_file.close() 2237 2238 # Check for empty files. 2239 if pdb_lines == []: 2240 raise RelaxError("The PDB file is empty.") 2241 2242 # Pre-process the lines, fixing PDB violations. 2243 pdb_lines = self._validate_records(pdb_lines) 2244 2245 # Secondary structure data. 2246 helices = [] 2247 sheets = [] 2248 2249 # Process the different sections. 2250 pdb_lines = self._parse_pdb_title(pdb_lines) 2251 pdb_lines = self._parse_pdb_prim_struct(pdb_lines) 2252 pdb_lines = self._parse_pdb_hetrogen(pdb_lines) 2253 pdb_lines = self._parse_pdb_ss(pdb_lines, read_mol=read_mol, helices=helices, sheets=sheets) 2254 pdb_lines = self._parse_pdb_connectivity_annotation(pdb_lines) 2255 pdb_lines = self._parse_pdb_misc(pdb_lines) 2256 pdb_lines = self._parse_pdb_transform(pdb_lines) 2257 2258 # Loop over all models in the PDB file. 2259 model_index = 0 2260 orig_model_num = [] 2261 mol_conts = [] 2262 orig_mol_num = [] 2263 for model_num, model_records in self._parse_pdb_coord(pdb_lines): 2264 # Only load the desired model. 2265 if read_model and model_num not in read_model: 2266 continue 2267 2268 # Store the original model number. 2269 orig_model_num.append(model_num) 2270 2271 # Loop over the molecules of the model. 2272 mol_conts.append([]) 2273 orig_mol_num.append([]) 2274 mol_index = 0 2275 new_mol_name = [] 2276 reset_serial = True 2277 for mol_num, mol_records in self._parse_mols_pdb(model_records): 2278 # Only load the desired model. 2279 if read_mol and mol_num not in read_mol: 2280 continue 2281 2282 # Set the target molecule name. 2283 if set_mol_name: 2284 # Single name. 2285 if len(set_mol_name) == 1: 2286 new_mol_name.append(set_mol_name[0]) 2287 2288 # Multiple names. 2289 else: 2290 new_mol_name.append(set_mol_name[mol_index]) 2291 2292 # Auto-generated molecule name. 2293 else: 2294 # Number of structures already present for the model. 2295 num_struct = 0 2296 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)): 2297 num_struct = len(self.structural_data[0].mol) 2298 2299 # Set the name to the file name plus the structure number. 2300 new_mol_name.append(file_root(file) + '_mol' + repr(mol_num+num_struct)) 2301 2302 # Generate the molecule container. 2303 mol = MolContainer() 2304 2305 # Fill the molecular data object. 2306 mol.fill_object_from_pdb(mol_records, alt_loc_select=alt_loc, reset_serial=reset_serial) 2307 2308 # Store the molecule container. 2309 mol_conts[model_index].append(mol) 2310 2311 # Store the original molecule number. 2312 orig_mol_num[model_index].append(mol_num) 2313 2314 # Increment the molecule index. 2315 mol_index = mol_index + 1 2316 2317 # No longer reset the serial number to 1. 2318 reset_serial = False 2319 2320 # Increment the model index. 2321 model_index = model_index + 1 2322 2323 # No data, so throw a warning and exit. 2324 if not len(mol_conts): 2325 warn(RelaxWarning("No structural data could be read from the file '%s'." % file_path)) 2326 return False 2327 2328 # Check the validity of the molecule names. 2329 invalid_name = [] 2330 for name in new_mol_name: 2331 if '#' in name: 2332 invalid_name.append(repr(name)) 2333 plural = '' 2334 if len(invalid_name) > 1: 2335 plural = 's' 2336 if len(invalid_name): 2337 raise RelaxError("The molecule name%s %s cannot contain the '#' character, as this is used in the spin ID strings as the molecule identifier." % (plural, human_readable_list(invalid_name))) 2338 2339 # Create the structural data data structures. 2340 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, helices=helices, sheets=sheets, verbosity=verbosity) 2341 2342 # Loading worked. 2343 return True

2344 2345

2346 - def load_xyz(self, file_path, read_mol=None, set_mol_name=None, read_model=None, set_model_num=None, verbosity=False):

2347 """Method for loading structures from a XYZ file. 2348 2349 @param file_path: The full path of the XYZ file. 2350 @type file_path: str 2351 @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. 2352 @type read_mol: None, int, or list of int 2353 @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. 2354 @type set_mol_name: None, str, or list of str 2355 @keyword read_model: The XYZ model to extract from the file. If set to None, then all models will be loaded. 2356 @type read_model: None, int, or list of int 2357 @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. 2358 @type set_model_num: None, int, or list of int 2359 @keyword verbosity: A flag which if True will cause messages to be printed. 2360 @type verbosity: bool 2361 @return: The status of the loading of the XYZ file. 2362 @rtype: bool 2363 """ 2364 2365 # Initial printout. 2366 if verbosity: 2367 print("\nInternal relax XYZ parser.") 2368 2369 # Test if the file exists. 2370 if not access(file_path, F_OK): 2371 # Exit indicating failure. 2372 return False 2373 2374 # Separate the file name and path. 2375 path, file = os.path.split(file_path) 2376 2377 # The absolute path. 2378 path_abs = abspath(curdir) + sep + path 2379 2380 # Convert the structure reading args into lists. 2381 if read_mol and not isinstance(read_mol, list): 2382 read_mol = [read_mol] 2383 if set_mol_name and not isinstance(set_mol_name, list): 2384 set_mol_name = [set_mol_name] 2385 if read_model and not isinstance(read_model, list): 2386 read_model = [read_model] 2387 if set_model_num and not isinstance(set_model_num, list): 2388 set_model_num = [set_model_num] 2389 2390 # Loop over all models in the XYZ file. 2391 mol_index = 0 2392 model_index = 0 2393 xyz_model_increment = 0 2394 orig_model_num = [] 2395 mol_conts = [] 2396 orig_mol_num = [] 2397 new_mol_name = [] 2398 for model_records in self._parse_models_xyz(file_path, verbosity=verbosity): 2399 # Increment the xyz_model_increment 2400 xyz_model_increment = xyz_model_increment +1 2401 2402 # Only load the desired model. 2403 if read_model and xyz_model_increment not in read_model: 2404 continue 2405 2406 # Store the original model number. 2407 orig_model_num.append(model_index) 2408 2409 # Loop over the molecules of the model. 2410 if read_mol and mol_index not in read_mol: 2411 continue 2412 2413 # Set the target molecule name. 2414 if set_mol_name: 2415 new_mol_name.append(set_mol_name[mol_index]) 2416 else: 2417 if mol_index == 0: 2418 #Set the name to the file name plus the structure number. 2419 new_mol_name.append(file_root(file) + '_mol' + repr(mol_index+1)) 2420 2421 # Store the original mol number. 2422 orig_mol_num.append(mol_index) 2423 2424 # Generate the molecule container. 2425 mol = MolContainer() 2426 2427 # Fill the molecular data object. 2428 mol.fill_object_from_xyz(model_records) 2429 2430 # Store the molecule container. 2431 mol_conts.append([]) 2432 mol_conts[model_index].append(mol) 2433 2434 # Increment the molecule index. 2435 mol_index = mol_index + 1 2436 2437 # Increment the model index. 2438 model_index = model_index + 1 2439 2440 # Check the validity of the molecule names. 2441 invalid_name = [] 2442 for name in new_mol_name: 2443 if '#' in name: 2444 invalid_name.append(repr(name)) 2445 plural = '' 2446 if len(invalid_name) > 1: 2447 plural = 's' 2448 if len(invalid_name): 2449 raise RelaxError("The molecule name%s %s cannot contain the '#' character, as this is used in the spin ID strings as the molecule identifier." % (plural, human_readable_list(invalid_name))) 2450 2451 # Create the structural data data structures. 2452 orig_mol_num = [[0]] 2453 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) 2454 2455 # Loading worked. 2456 return True

2457 2458

2459 - def mean(self):

2460 """Calculate the mean structure from all models in the structural data object.""" 2461 2462 # Create a new model for the mean structure. 2463 num = self.num_models() 2464 self.add_model() 2465 mean_model = self.structural_data[-1] 2466 2467 # The selection object. 2468 selection = self.selection() 2469 2470 # Loop over the molecules and atoms. 2471 for mol_index, i in selection.loop(): 2472 # Set the mean structure coordinate to zero. 2473 mean_model.mol[mol_index].x[i] = 0.0 2474 mean_model.mol[mol_index].y[i] = 0.0 2475 mean_model.mol[mol_index].z[i] = 0.0 2476 2477 # Loop over the models and sum the coordinates. 2478 for model_index in range(num): 2479 model_cont = self.structural_data[model_index] 2480 mean_model.mol[mol_index].x[i] += model_cont.mol[mol_index].x[i] 2481 mean_model.mol[mol_index].y[i] += model_cont.mol[mol_index].y[i] 2482 mean_model.mol[mol_index].z[i] += model_cont.mol[mol_index].z[i] 2483 2484 # Averages. 2485 mean_model.mol[mol_index].x[i] /= num 2486 mean_model.mol[mol_index].y[i] /= num 2487 mean_model.mol[mol_index].z[i] /= num 2488 2489 # Delete all models but the mean. 2490 for model_index in reversed(list(range(num))): 2491 self.delete(model=self.structural_data[model_index].num)

2492 2493

2494 - def model_list(self):

2495 """Create a list of all models. 2496 2497 @return: The list of all models. 2498 @rtype: list of int 2499 """ 2500 2501 # Assemble the list. 2502 models = [] 2503 for model in self.model_loop(): 2504 models.append(model.num) 2505 2506 # Return the list. 2507 return models

2508 2509

2510 - def model_loop(self, model=None):

2511 """Generator method for looping over the models in numerical order. 2512 2513 @keyword model: Limit the loop to a single number. 2514 @type model: int 2515 @return: The model structural object. 2516 @rtype: ModelContainer container 2517 """ 2518 2519 # A single model. 2520 if model: 2521 for i in range(len(self.structural_data)): 2522 if self.structural_data[i].num == model: 2523 yield self.structural_data[i] 2524 2525 # All models. 2526 else: 2527 # The models. 2528 model_nums = [] 2529 for i in range(len(self.structural_data)): 2530 if self.structural_data[i].num != None: 2531 model_nums.append(self.structural_data[i].num) 2532 2533 # Sort. 2534 if model_nums: 2535 model_nums.sort() 2536 2537 # Loop over the models in order. 2538 for model_num in model_nums: 2539 # Find the model. 2540 for i in range(len(self.structural_data)): 2541 # Yield the model. 2542 if self.structural_data[i].num == model_num: 2543 yield self.structural_data[i] 2544 2545 # No models, so just yield the single container. 2546 if not model_nums: 2547 yield self.structural_data[0]

2548 2549

2550 - def num_models(self):

2551 """Method for returning the number of models. 2552 2553 @return: The number of models in the structural object. 2554 @rtype: int 2555 """ 2556 2557 return len(self.structural_data)

2558 2559

2560 - def num_molecules(self):

2561 """Method for returning the number of molecules. 2562 2563 @return: The number of molecules in the structural object. 2564 @rtype: int 2565 """ 2566 2567 # No data. 2568 if self.empty(): 2569 return 0 2570 2571 # Validate the structural object. 2572 self.validate() 2573 2574 # Return the number. 2575 return len(self.structural_data[0].mol)

2576 2577

2578 - def one_letter_codes(self, mol_name=None, selection=None):

2579 """Generate and return the one letter code sequence for the given molecule. 2580 2581 @keyword mol_name: The name of the molecule to return the one letter codes for. 2582 @type mol_name: str 2583 @keyword selection: The internal structural selection object. This is obtained by calling the selection() method with the atom ID string. 2584 @type selection: lib.structure.internal.Internal_selection instance 2585 @return: The one letter code sequence for the given molecule. 2586 @rtype: str 2587 """ 2588 2589 # Initialise. 2590 codes = '' 2591 2592 # Use the first model. 2593 model = self.structural_data[0] 2594 2595 # Residue numbers. 2596 res_nums = [] 2597 2598 # Loop over the molecules. 2599 for mol_index, i in selection.loop(): 2600 # Alias. 2601 mol = model.mol[mol_index] 2602 2603 # Skip non-matching molecules. 2604 if mol_name and mol_name != mol.mol_name: 2605 continue 2606 2607 # Not a new residue. 2608 if mol.res_num[i] in res_nums: 2609 continue 2610 2611 # Convert to the one letter code and store the residue number. 2612 codes += aa_codes_three_to_one(mol.res_name[i]) 2613 res_nums.append(mol.res_num[i]) 2614 2615 # Return the codes. 2616 return codes

2617 2618

2619 - 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, helices=None, sheets=None, verbosity=1, merge=False):

2620 """From the given structural data, expand the structural data data structure. 2621 2622 @param data_matrix: A matrix of structural objects. The first dimension is the models, the second is the molecules. 2623 @type data_matrix: list of lists of structural objects 2624 @keyword orig_model_num: The original model numbers (for storage). 2625 @type orig_model_num: list of int 2626 @keyword set_model_num: The new model numbers (for model renumbering). 2627 @type set_model_num: list of int 2628 @keyword orig_mol_num: The original molecule numbers (for storage). The dimensions match the data matrix. 2629 @type orig_mol_num: list of list of int 2630 @keyword set_mol_name: The new molecule names. 2631 @type set_mol_name: list of str 2632 @keyword file_name: The name of the file from which the molecular data has been extracted. 2633 @type file_name: None or str 2634 @keyword file_path: The full path to the file specified by 'file_name'. 2635 @type file_path: None or str 2636 @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. 2637 @type file_path_abs: None or str 2638 @keyword helices: The helix secondary structure information to store. The first dimension corresponds to each helix element. The second dimension consists of the helix_id, mol_init_index, init_res_name, init_seq_num, mol_end_index, end_res_name, end_seq_num, helix_class, and length. 2639 @type helices: list of lists 2640 @keyword sheets: The sheet secondary structure information to store. 2641 @type sheets: list of lists 2642 @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. 2643 @type verbosity: int 2644 @keyword merge: A flag which if set to True will try to merge the structure into the currently loaded structures. 2645 @type merge: bool 2646 """ 2647 2648 # Test the number of models. 2649 if len(orig_model_num) != len(data_matrix): 2650 raise RelaxError("Structural data mismatch, %s original models verses %s in the structural data." % (len(orig_model_num), len(data_matrix))) 2651 2652 # Test the number of molecules. 2653 for i in range(len(data_matrix)): 2654 if len(orig_mol_num[i]) != len(data_matrix[i]): 2655 raise RelaxError("Structural data mismatch, %s original molecules verses %s in the structural data." % (len(orig_mol_num[i]), len(data_matrix[i]))) 2656 2657 # Model numbers do not change. 2658 if not set_model_num: 2659 set_model_num = orig_model_num 2660 2661 # Test the model mapping. 2662 if len(set_model_num) != len(data_matrix): 2663 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))) 2664 2665 # Test the molecule mapping. 2666 if len(set_mol_name) != len(data_matrix[0]): 2667 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]))) 2668 2669 # Test that the target models and structures are absent, and get the already existing model numbers. 2670 for i in range(len(set_model_num)): 2671 # Merging flag is set, so skip the checks. 2672 if merge: 2673 continue 2674 2675 # A new model, so no need to check. 2676 if not set_model_num[i] in self.structural_data.current_models: 2677 continue 2678 2679 # Loop over the structures. 2680 index = self.structural_data.current_models.index(set_model_num[i]) 2681 for j in range(len(self.structural_data[index].mol)): 2682 if self.structural_data[index].mol[j].mol_name in set_mol_name: 2683 raise RelaxError("The molecule '%s' of model %s already exists." % (self.structural_data[i].mol[j].mol_name, self.structural_data[i].num)) 2684 2685 # Loop over the models. 2686 store_mol_num_orig = [] 2687 store_mol_num_new = [] 2688 for i in range(len(set_model_num)): 2689 store_mol_num_orig.append([]) 2690 store_mol_num_new.append([]) 2691 2692 # The model doesn't currently exist. 2693 if set_model_num[i] not in self.structural_data.current_models: 2694 # Create the model. 2695 self.structural_data.add_item(set_model_num[i]) 2696 2697 # Get the model. 2698 model = self.structural_data[-1] 2699 2700 # Otherwise get the pre-existing model. 2701 else: 2702 model = self.structural_data[self.structural_data.current_models.index(set_model_num[i])] 2703 2704 # Loop over the molecules. 2705 for j in range(len(set_mol_name)): 2706 # Override the merge argument if the molecule does not exist. 2707 found = False 2708 merge_new = True 2709 for k in range(len(model.mol)): 2710 if model.mol[k].mol_name == set_mol_name[j]: 2711 found = True 2712 index = k 2713 if not found: 2714 merge_new = False 2715 2716 # Printout. 2717 if verbosity: 2718 # Model text formatting. 2719 orig_model_text = '' 2720 if orig_model_num[i] != None: 2721 orig_model_text = " of model %s" % orig_model_num[i] 2722 new_model_text = '' 2723 if set_model_num[i] != None: 2724 if merge_new: 2725 new_model_text += ' of' 2726 else: 2727 new_model_text += ' to' 2728 new_model_text += ' model %s' % set_model_num[i] 2729 2730 # The full text. 2731 if merge_new: 2732 print("Merging with molecule '%s'%s (from the original molecule number %s%s)." % (set_mol_name[j], new_model_text, orig_mol_num[i][j], orig_model_text)) 2733 else: 2734 print("Adding molecule '%s'%s (from the original molecule number %s%s)." % (set_mol_name[j], new_model_text, orig_mol_num[i][j], orig_model_text)) 2735 2736 # The index of the new molecule to add or merge. 2737 if not merge_new: 2738 index = len(model.mol) 2739 2740 # Store the index+1 as the new molecule number, and store the original molecule number. 2741 store_mol_num_new[i].append(index+1) 2742 store_mol_num_orig[i].append(orig_mol_num[i][j]) 2743 2744 # Consistency check. 2745 if model.num != self.structural_data[0].num and self.structural_data[0].mol[index].mol_name != set_mol_name[j]: 2746 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)) 2747 2748 # Pack the structures. 2749 if merge_new: 2750 mol = model.mol.merge_item(mol_name=set_mol_name[j], mol_cont=data_matrix[i][j]) 2751 else: 2752 mol = model.mol.add_item(mol_name=set_mol_name[j], mol_cont=data_matrix[i][j]) 2753 2754 # Set the molecule name and store the structure file info. 2755 mol.mol_name = set_mol_name[j] 2756 mol.file_name = file_name 2757 mol.file_path = file_path 2758 mol.file_path_abs = file_path_abs 2759 mol.file_mol_num = orig_mol_num[i][j] 2760 mol.file_model = orig_model_num[i] 2761 2762 # Sort the structural data if a merge occurred. 2763 mol._sort() 2764 2765 # Consistency checks. 2766 for model_index in range(len(store_mol_num_new)): 2767 if len(store_mol_num_new[model_index]) != len(store_mol_num_orig[model_index]): 2768 raise RelaxFault 2769 2770 # The helix secondary structure information. 2771 if helices != None: 2772 # Initialise the data structure if required. 2773 if not hasattr(self, 'helices'): 2774 self.helices = [] 2775 2776 # Loop over each helix. 2777 for i in range(len(helices)): 2778 # Append the data. 2779 self.helices.append(helices[i]) 2780 2781 # The initial molecule index remapping. 2782 mol_init_index = store_mol_num_new[0][store_mol_num_orig[0].index(helices[i][1]+1)] - 1 2783 self.helices[-1][1] = mol_init_index 2784 2785 # The end molecule index remapping. 2786 mol_end_index = store_mol_num_new[0][store_mol_num_orig[0].index(helices[i][4]+1)] - 1 2787 self.helices[-1][4] = mol_end_index 2788 2789 # The sheet secondary structure information. 2790 if sheets != None: 2791 # Initialise the data structure if required. 2792 if not hasattr(self, 'sheets'): 2793 self.sheets = [] 2794 2795 # Loop over each sheet. 2796 for i in range(len(sheets)): 2797 # Append the data. 2798 self.sheets.append(sheets[i]) 2799 2800 # The initial molecule index remappings. 2801 mol_init_index = store_mol_num_new[0][store_mol_num_orig[0].index(sheets[i][4]+1)] - 1 2802 self.sheets[-1][4] = mol_init_index 2803 2804 2805 # The end molecule index remapping. 2806 mol_end_index = store_mol_num_new[0][store_mol_num_orig[0].index(sheets[i][8]+1)] - 1 2807 self.sheets[-1][8] = mol_end_index 2808 2809 # The current molecule index remapping. 2810 if sheets[i][14] != None: 2811 mol_cur_index = store_mol_num_new[0][store_mol_num_orig[0].index(sheets[i][14]+1)] - 1 2812 self.sheets[-1][14] = mol_cur_index 2813 2814 # The previous molecule index remapping. 2815 if sheets[i][19] != None: 2816 mol_prev_index = store_mol_num_new[0][store_mol_num_orig[0].index(sheets[i][19]+1)] - 1 2817 self.sheets[-1][19] = mol_prev_index

2818 2819

2820 - def rotate(self, R=None, origin=None, model=None, selection=None):

2821 """Rotate the structural information about the given origin. 2822 2823 @keyword R: The forwards rotation matrix. 2824 @type R: numpy 3D, rank-2 array 2825 @keyword origin: The origin of the rotation. 2826 @type origin: numpy 3D, rank-1 array 2827 @keyword model: The model to rotate. If None, all models will be rotated. 2828 @type model: int 2829 @keyword selection: The internal structural selection object. This is obtained by calling the selection() method with the atom ID string. 2830 @type selection: lib.structure.internal.Internal_selection instance 2831 """ 2832 2833 # Loop over the models. 2834 for model_cont in self.model_loop(model): 2835 # Loop over all molecules and atoms in the selection. 2836 for mol_index, i in selection.loop(): 2837 mol = model_cont.mol[mol_index] 2838 2839 # The origin to atom vector. 2840 vect = array([mol.x[i], mol.y[i], mol.z[i]], float64) - origin 2841 2842 # Rotation. 2843 rot_vect = dot(R, vect) 2844 2845 # The new position. 2846 pos = rot_vect + origin 2847 mol.x[i] = pos[0] 2848 mol.y[i] = pos[1] 2849 mol.z[i] = pos[2]

2850 2851

2852 - def selection(self, atom_id=None, inv=False):

2853 """Convert the atom ID string into a special internal selection object for speed. 2854 2855 @keyword atom_id: The molecule, residue, and atom identifier string. Only atoms matching this selection will be used. 2856 @type atom_id: str or None 2857 @keyword inv: A flag which if True will cause the selection to be inverted. 2858 @type inv: bool 2859 @return: The internal structural selection object. 2860 @rtype: Internal_selection instance 2861 """ 2862 2863 # Initialise the internal structural selection object. 2864 selection = Internal_selection() 2865 2866 # Split up the atom ID, to see if a molecule has been specified. 2867 mol_token, res_token, spin_token = tokenise(atom_id) 2868 2869 # Generate the atom ID selection object. 2870 sel_obj = None 2871 if atom_id: 2872 sel_obj = Selection(atom_id) 2873 2874 # Obtain all data from the first model (except the position data). 2875 model = self.structural_data[0] 2876 2877 # Loop over the molecules. 2878 for mol_index in range(len(model.mol)): 2879 mol = model.mol[mol_index] 2880 2881 # Skip non-matching molecules. 2882 if not inv: 2883 if sel_obj and not sel_obj.contains_mol(mol.mol_name): 2884 continue 2885 2886 # Skip matching molecules. 2887 else: 2888 if (not sel_obj) or (mol_token != None and sel_obj.contains_mol(mol.mol_name)): 2889 continue 2890 2891 # Add the molecule index. 2892 selection.add_mol(mol_index=mol_index) 2893 2894 # Loop over the atoms. 2895 for i in range(len(mol.atom_num)): 2896 # Skip non-matching atoms. 2897 if not inv: 2898 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): 2899 continue 2900 2901 # Skip matching atoms. 2902 else: 2903 if (not sel_obj) or sel_obj.contains_spin(mol.atom_num[i], mol.atom_name[i], mol.res_num[i], mol.res_name[i], mol.mol_name): 2904 continue 2905 2906 # Add the atom index. 2907 selection.add_atom(mol_index=mol_index, atom_index=i) 2908 2909 # Return the object. 2910 return selection

2911 2912

2913 - def set_model(self, model_orig=None, model_new=None):

2914 """Set or reset the model number. 2915 @keyword model_orig: The original model number. Leave as None if no models are currently present. 2916 @type model_orig: None or int 2917 @keyword model_new: The new model number to set the model to. 2918 @type model_new: int 2919 """ 2920 2921 # Check. 2922 if model_orig == None and self.num_models() != 1: 2923 raise RelaxError("If the original model number is not supplied, only one model in the current structural object is allowed, but %s were found." % self.num_models()) 2924 2925 # Set the single model number. 2926 if model_orig == None: 2927 self.structural_data[0].num = model_new 2928 self.structural_data.current_models[0] = model_new 2929 return 2930 2931 # Find the model and set the number. 2932 set = False 2933 for i in range(len(self.structural_data)): 2934 if model_orig == self.structural_data[i].num: 2935 self.structural_data[i].num = model_new 2936 self.structural_data.current_models[i] = model_new 2937 set = True 2938 2939 # Sanity check. 2940 if not set: 2941 raise RelaxError("The original model number %s could not be found in the structural object." % model_orig)

2942 2943

2944 - def target_mol_name(self, set=None, target=None, index=None, mol_num=None, file=None):

2945 """Add the new molecule name to the target data structure. 2946 2947 @keyword set: The list of new molecule names. If not supplied, the names will be generated from the file name. 2948 @type set: None or list of str 2949 @keyword target: The target molecule name data structure to which the new name will be appended. 2950 @type target: list 2951 @keyword index: The molecule index, matching the set argument. 2952 @type index: int 2953 @keyword mol_num: The molecule number. 2954 @type mol_num: int 2955 @keyword file: The name of the file, excluding all directories. 2956 @type file: str 2957 """ 2958 2959 # Set the target molecule name. 2960 if set: 2961 target.append(set[index]) 2962 else: 2963 # Set the name to the file name plus the structure number. 2964 target.append(file_root(file) + '_mol' + repr(mol_num))

2965 2966

2967 - def translate(self, T=None, model=None, selection=None):

2968 """Displace the structural information by the given translation vector. 2969 2970 @keyword T: The translation vector. 2971 @type T: numpy 3D, rank-1 array 2972 @keyword model: The model to rotate. If None, all models will be rotated. 2973 @type model: int 2974 @keyword selection: The internal structural selection object. This is obtained by calling the selection() method with the atom ID string. 2975 @type selection: lib.structure.internal.Internal_selection instance 2976 """ 2977 2978 # Loop over the models. 2979 for model_cont in self.model_loop(model): 2980 # Loop over all molecules and atoms in the selection. 2981 for mol_index, i in selection.loop(): 2982 mol = model_cont.mol[mol_index] 2983 2984 # Translate. 2985 mol.x[i] = mol.x[i] + T[0] 2986 mol.y[i] = mol.y[i] + T[1] 2987 mol.z[i] = mol.z[i] + T[2]

2988 2989

2990 - def to_xml(self, doc, element):

2991 """Prototype method for converting the structural object to an XML representation. 2992 2993 @param doc: The XML document object. 2994 @type doc: xml.dom.minidom.Document instance 2995 @param element: The element to add the alignment tensors XML element to. 2996 @type element: XML element object 2997 """ 2998 2999 # Create the structural element and add it to the higher level element. 3000 str_element = doc.createElement('structure') 3001 element.appendChild(str_element) 3002 3003 # Set the structural attributes. 3004 str_element.setAttribute('desc', 'Structural information') 3005 3006 # No contents to store, so pack up the structural containers. 3007 if not self.structural_data.is_empty(): 3008 self.structural_data.to_xml(doc, str_element) 3009 3010 # The structural metadata. 3011 metadata = ['helices', 'sheets'] 3012 for name in metadata: 3013 # The metadata does not exist. 3014 if not hasattr(self, name): 3015 continue 3016 3017 # Get the object. 3018 obj = getattr(self, name) 3019 3020 # Create a new element for this object, and add it to the main element. 3021 sub_elem = doc.createElement(name) 3022 str_element.appendChild(sub_elem) 3023 3024 # Add the value to the sub element. 3025 object_to_xml(doc, sub_elem, value=obj) 3026 3027 # The displacement structure. 3028 if hasattr(self, 'displacements'): 3029 # Create an XML element. 3030 disp_element = doc.createElement('displacements') 3031 str_element.appendChild(disp_element) 3032 3033 # Set the attributes. 3034 disp_element.setAttribute('desc', 'The rotational and translational displacements between models') 3035 3036 # Add the displacement data. 3037 self.displacements.to_xml(doc, disp_element) 3038 3039 # Blacklisted objects. 3040 blacklist = ['structural_data', 'displacements'] + metadata + list(Internal.__dict__.keys()) + list(self.__class__.__dict__.keys()) + list(object.__dict__.keys()) 3041 3042 # Add all simple python objects within the container to the XML element. 3043 fill_object_contents(doc, str_element, object=self, blacklist=blacklist)

3044 3045

3046 - def validate_models(self, verbosity=1):

3047 """Check that the models are consistent with each other. 3048 3049 This checks that the primary structure is identical between the models. 3050 3051 3052 @keyword verbosity: If 0, then all printouts will be silenced. 3053 @type verbosity: int 3054 """ 3055 3056 # Print out. 3057 if verbosity: 3058 print("Validating models:") 3059 3060 # Loop over the models. 3061 for i in range(len(self.structural_data)): 3062 # Check the molecules. 3063 if len(self.structural_data[0].mol) != len(self.structural_data[i].mol): 3064 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))) 3065 3066 # Loop over the molecules. 3067 for j in range(len(self.structural_data[i].mol)): 3068 # Alias the molecules. 3069 mol = self.structural_data[i].mol[j] 3070 mol_ref = self.structural_data[0].mol[j] 3071 3072 # Check the names. 3073 if mol.mol_name != mol_ref.mol_name: 3074 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)) 3075 3076 # Loop over the atoms. 3077 for k in range(len(mol.atom_name)): 3078 # Check all data. 3079 same = True 3080 if mol.chain_id[k] != mol_ref.chain_id[k]: 3081 same = False 3082 elif mol.atom_num[k] != mol_ref.atom_num[k]: 3083 same = False 3084 elif mol.atom_name[k] != mol_ref.atom_name[k]: 3085 same = False 3086 elif mol.res_name[k] != mol_ref.res_name[k]: 3087 same = False 3088 elif mol.res_num[k] != mol_ref.res_num[k]: 3089 same = False 3090 elif mol.seg_id[k] != mol_ref.seg_id[k]: 3091 same = False 3092 elif mol.element[k] != mol_ref.element[k]: 3093 same = False 3094 3095 # Failure. 3096 if not same: 3097 print("%-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]), '')) 3098 print("%-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]), '')) 3099 raise RelaxError("The atoms of model %s do not match the first model." % self.structural_data[i].num) 3100 3101 # Final printout. 3102 if verbosity: 3103 print("\tAll models are consistent")

3104 3105

3106 - def write_pdb(self, file, model_num=None):

3107 """Method for the creation of a PDB file from the structural data. 3108 3109 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: 3110 3111 0. Residue number. 3112 1. Residue name. 3113 2. Chain ID. 3114 3. Total number of atoms in the residue. 3115 4. Number of H atoms in the residue. 3116 5. Number of C atoms in the residue. 3117 3118 3119 @param file: The PDB file object. This object must be writable. 3120 @type file: file object 3121 @keyword model_num: The model to place into the PDB file. If not supplied, then all models will be placed into the file. 3122 @type model_num: None or int 3123 """ 3124 3125 # Validate the structural data. 3126 self.validate() 3127 3128 # Print out. 3129 print("\nCreating the PDB records\n") 3130 3131 # Write some initial remarks. 3132 print("REMARK") 3133 pdb_write.remark(file, num=4, remark="This file complies with format v. 3.30, Jul-2011.") 3134 pdb_write.remark(file, num=40, remark=None) 3135 pdb_write.remark(file, num=40, remark="Created using relax (http://www.nmr-relax.com).") 3136 pdb_write.remark(file, num=40, remark=None) 3137 if RELAX_VERSION: 3138 pdb_write.remark(file, num=40, remark="relax version %s." % RELAX_VERSION) 3139 pdb_write.remark(file, num=40, remark="Created on %s." % asctime()) 3140 num_remark = 2 3141 3142 # Determine if model records will be created. 3143 model_records = False 3144 for model in self.model_loop(): 3145 if hasattr(model, 'num') and model.num != None: 3146 model_records = True 3147 3148 3149 #################### 3150 # Hetrogen section # 3151 #################### 3152 3153 # Initialise the hetrogen info array. 3154 het_data = [] 3155 het_data_coll = [] 3156 3157 # Loop over the molecules of the first model. 3158 index = 0 3159 for mol in self.structural_data[0].mol_loop(): 3160 # Check the validity of the data. 3161 self._validate_data_arrays(mol) 3162 3163 # Append an empty array for this molecule. 3164 het_data.append([]) 3165 3166 # Collect the non-standard residue info. 3167 for i in range(len(mol.atom_name)): 3168 # Skip non-HETATM records and HETATM records with no residue info. 3169 if mol.pdb_record[i] != 'HETATM' or mol.res_name[i] == None: 3170 continue 3171 3172 # Skip waters. 3173 if mol.res_name[i] == 'HOH': 3174 continue 3175 3176 # If the residue is not already stored initialise a new het_data element. 3177 # (residue number, residue name, chain ID, number of atoms, atom count array). 3178 if not het_data[index] or not mol.res_num[i] == het_data[index][-1][0]: 3179 het_data[index].append([mol.res_num[i], mol.res_name[i], CHAIN_ID_LIST[index], 0, []]) 3180 3181 # Catch missing chain_ids. 3182 if het_data[index][-1][2] == None: 3183 het_data[index][-1][2] = '' 3184 3185 # Total atom count. 3186 het_data[index][-1][3] = het_data[index][-1][3] + 1 3187 3188 # Find if the atom has already a count entry. 3189 entry = False 3190 for j in range(len(het_data[index][-1][4])): 3191 if mol.element[i] == het_data[index][-1][4][j][0]: 3192 entry = True 3193 3194 # Create a new specific atom count entry. 3195 if not entry: 3196 het_data[index][-1][4].append([mol.element[i], 0]) 3197 3198 # Increment the specific atom count. 3199 for j in range(len(het_data[index][-1][4])): 3200 if mol.element[i] == het_data[index][-1][4][j][0]: 3201 het_data[index][-1][4][j][1] = het_data[index][-1][4][j][1] + 1 3202 3203 # Create the collective hetrogen info data structure. 3204 for i in range(len(het_data[index])): 3205 # Find the entry in the collective structure. 3206 found = False 3207 for j in range(len(het_data_coll)): 3208 # Different residue numbers. 3209 if het_data[index][i][0] != het_data_coll[j][0]: 3210 continue 3211 3212 # Different chain IDs. 3213 if het_data_coll[j][2] != het_data[index][i][2]: 3214 continue 3215 3216 # Change the flag. 3217 found = True 3218 3219 # If there is no match, add the new residue to the collective. 3220 if not found: 3221 het_data_coll.append(het_data[index][i]) 3222 3223 # Increment the molecule index. 3224 index = index + 1 3225 3226 3227 # The HET records. 3228 ################## 3229 3230 # Print out. 3231 print("HET") 3232 3233 # Write the HET records. 3234 for het in het_data_coll: 3235 pdb_write.het(file, het_id=het[1], chain_id=het[2], seq_num=het[0], num_het_atoms=het[3]) 3236 3237 3238 # The HETNAM records. 3239 ##################### 3240 3241 # Print out. 3242 print("HETNAM") 3243 3244 # Loop over the non-standard residues. 3245 residues = [] 3246 for het in het_data_coll: 3247 # Test if the residue HETNAM record as already been written (otherwise store its name). 3248 if het[1] in residues: 3249 continue 3250 else: 3251 residues.append(het[1]) 3252 3253 # Get the chemical name. 3254 chemical_name = self._get_chemical_name(het[1]) 3255 if not chemical_name: 3256 chemical_name = 'Unknown' 3257 3258 # Write the HETNAM records. 3259 pdb_write.hetnam(file, het_id=het[1], text=chemical_name) 3260 3261 3262 # The FORMUL records. 3263 ##################### 3264 3265 # Print out. 3266 print("FORMUL") 3267 3268 # Loop over the non-standard residues and generate and write the chemical formula. 3269 residues = [] 3270 for i in range(len(het_data_coll)): 3271 # Alias. 3272 het = het_data_coll[i] 3273 3274 # Test if the residue HETNAM record as already been written (otherwise store its name). 3275 if het[1] in residues: 3276 continue 3277 else: 3278 residues.append(het[1]) 3279 3280 # Initialise the chemical formula. 3281 formula = '' 3282 3283 # Loop over the atoms. 3284 for atom_count in het[4]: 3285 formula = formula + atom_count[0] + repr(atom_count[1]) 3286 3287 # The FORMUL record (chemical formula). 3288 pdb_write.formul(file, comp_num=i+1, het_id=het[1], text=formula) 3289 3290 3291 ############################### 3292 # Secondary structure section # 3293 ############################### 3294 3295 # The HELIX records. 3296 #################### 3297 3298 if hasattr(self, 'helices') and len(self.helices): 3299 # Printout. 3300 print("HELIX") 3301 3302 # Loop over and unpack the helix data. 3303 index = 1 3304 for helix_id, mol_init_index, init_res_name, init_seq_num, mol_end_index, end_res_name, end_seq_num, helix_class, length in self.helices: 3305 # The chain IDs. 3306 if mol_init_index == None: 3307 mol_init_index = 0 3308 if mol_end_index == None: 3309 mol_end_index = 0 3310 3311 pdb_write.helix(file, ser_num=index, helix_id=helix_id, init_chain_id=CHAIN_ID_LIST[mol_init_index], init_res_name=init_res_name, init_seq_num=init_seq_num, end_chain_id=CHAIN_ID_LIST[mol_end_index], end_res_name=end_res_name, end_seq_num=end_seq_num, helix_class=helix_class, length=length) 3312 index += 1 3313 3314 # The SHEET records. 3315 #################### 3316 3317 if hasattr(self, 'sheets') and len(self.sheets): 3318 # Printout. 3319 print("SHEET") 3320 3321 # Loop over and unpack the helix data. 3322 index = 1 3323 for strand, sheet_id, num_strands, init_res_name, mol_init_index, init_seq_num, init_icode, end_res_name, mol_end_index, end_seq_num, end_icode, sense, cur_atom, cur_res_name, mol_cur_index, cur_res_seq, cur_icode, prev_atom, prev_res_name, mol_prev_index, prev_res_seq, prev_icode in self.sheets: 3324 # Translate molecule indices to chain IDs. 3325 if mol_init_index == None: 3326 mol_init_index = 0 3327 if mol_end_index == None: 3328 mol_end_index = 0 3329 init_chain_id = CHAIN_ID_LIST[mol_init_index] 3330 end_chain_id = CHAIN_ID_LIST[mol_end_index] 3331 cur_chain_id = None 3332 if mol_cur_index != None: 3333 cur_chain_id = CHAIN_ID_LIST[mol_cur_index] 3334 prev_chain_id = None 3335 if mol_prev_index != None: 3336 prev_chain_id = CHAIN_ID_LIST[mol_prev_index] 3337 3338 # Write out. 3339 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) 3340 index += 1 3341 3342 3343 ###################### 3344 # Coordinate section # 3345 ###################### 3346 3347 # Initial printout if models are present. 3348 if model_records: 3349 print("\nMODEL records:") 3350 3351 # Loop over the models. 3352 for model in self.model_loop(model_num): 3353 # Initialise record counts. 3354 num_hetatm = 0 3355 num_atom = 0 3356 num_ter = 0 3357 ser_num = 1 3358 3359 # MODEL record, for multiple models. 3360 #################################### 3361 3362 if model_records: 3363 # Printout. 3364 sys.stdout.write('.') 3365 3366 # Write the model record. 3367 pdb_write.model(file, serial=model.num) 3368 3369 3370 # Add the atomic coordinate records (ATOM, HETATM, and TER). 3371 ############################################################ 3372 3373 # Loop over the molecules. 3374 index = 0 3375 for mol in model.mol_loop(): 3376 # Printout. 3377 if not model_records: 3378 print("ATOM, HETATM, TER") 3379 3380 # Loop over the atomic data. 3381 atom_record = False 3382 for i in range(len(mol.atom_name)): 3383 # Write the ATOM record. 3384 if mol.pdb_record[i] in [None, 'ATOM']: 3385 atom_record = True 3386 3387 # Write out. 3388 pdb_write.atom(file, serial=ser_num, name=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, element=mol.element[i]) 3389 num_atom += 1 3390 ser_num += 1 3391 3392 # Info for the TER record. 3393 ter_num = ser_num + 1 3394 ter_name = mol.res_name[i] 3395 ter_chain_id = CHAIN_ID_LIST[index] 3396 ter_res_num = mol.res_num[i] 3397 3398 # Finish the ATOM section with the TER record. 3399 if atom_record: 3400 pdb_write.ter(file, serial=ser_num, res_name=ter_name, chain_id=CHAIN_ID_LIST[index], res_seq=ter_res_num) 3401 num_ter += 1 3402 ser_num += 1 3403 3404 # Loop over the atomic data. 3405 count_shift = False 3406 for i in range(len(mol.atom_name)): 3407 # Write the HETATM record. 3408 if mol.pdb_record[i] == 'HETATM': 3409 # Increment the atom number if a TER record was created. 3410 if atom_record and ser_num == ter_num: 3411 count_shift = True 3412 3413 # Write out. 3414 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]) 3415 num_hetatm += 1 3416 ser_num += 1 3417 3418 # Increment the molecule index. 3419 index += 1 3420 3421 3422 # ENDMDL record, for multiple structures. 3423 ######################################## 3424 3425 if model_records: 3426 if not model_records: 3427 print("ENDMDL") 3428 pdb_write.endmdl(file) 3429 3430 3431 # Create the CONECT records. 3432 ############################ 3433 3434 # Print out. 3435 if model_records: 3436 sys.stdout.write('\n') 3437 print("CONECT") 3438 3439 # Initialise record counts. 3440 num_conect = 0 3441 3442 # The per molecule incremented atom counts. 3443 atom_counts = [0] 3444 index = 0 3445 for mol in self.structural_data[0].mol_loop(): 3446 if index == 0: 3447 atom_counts.append(len(mol.atom_name)) 3448 else: 3449 atom_counts.append(atom_counts[index] + len(mol.atom_name)) 3450 index += 1 3451 3452 # Loop over the molecules of the first model. 3453 index = 0 3454 for mol in self.structural_data[0].mol_loop(): 3455 # Loop over the atoms. 3456 for i in range(len(mol.atom_name)): 3457 # No bonded atoms, hence no CONECT record is required. 3458 if not len(mol.bonded[i]): 3459 continue 3460 3461 # Initialise some data structures. 3462 flush = 0 3463 bonded_index = 0 3464 bonded = ['', '', '', ''] 3465 bonded_shifted = ['', '', '', ''] 3466 3467 # Loop over the bonded atoms. 3468 for j in range(len(mol.bonded[i])): 3469 # End of the array, hence create the CONECT record in this iteration. 3470 if j == len(mol.bonded[i])-1: 3471 flush = True 3472 3473 # Only four covalently bonded atoms allowed in one CONECT record. 3474 if bonded_index == 3: 3475 flush = True 3476 3477 # Get the bonded atom index. 3478 bonded[bonded_index] = mol.bonded[i][j] 3479 3480 # Increment the bonded_index value. 3481 bonded_index = bonded_index + 1 3482 3483 # Generate the CONECT record and increment the counter. 3484 if flush: 3485 # Convert the atom indices to atom numbers. 3486 for k in range(4): 3487 if bonded[k] != '': 3488 bonded_shifted[k] = bonded[k] + 1 + atom_counts[index] 3489 3490 # Write the CONECT record. 3491 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]) 3492 3493 # Reset the flush flag, the bonded atom count, and the bonded atom names. 3494 flush = False 3495 bonded_index = 0 3496 bonded = ['', '', '', ''] 3497 3498 # Increment the CONECT record count. 3499 num_conect = num_conect + 1 3500 3501 # Increment the molecule index. 3502 index += 1 3503 3504 3505 # MASTER record. 3506 ################ 3507 3508 print("\nMASTER") 3509 pdb_write.master(file, num_het=len(het_data_coll), num_coord=num_atom+num_hetatm, num_ter=num_ter, num_conect=num_conect) 3510 3511 3512 # END. 3513 ###### 3514 3515 print("END") 3516 pdb_write.end(file)

3517 3518

3519 - def validate(self):

3520 """Check the integrity of the structural data. 3521 3522 The number of molecules must be the same in all models. 3523 """ 3524 3525 # Reference number of molecules. 3526 num_mols = len(self.structural_data[0].mol) 3527 3528 # Loop over all other models. 3529 for i in range(1, len(self.structural_data)): 3530 # Model alias. 3531 model_i = self.structural_data[i] 3532 3533 # Size check. 3534 if num_mols != len(model_i.mol): 3535 raise RelaxError("The structural object is not valid - the number of molecules is not the same for all models.") 3536 3537 # Molecule name check. 3538 for j in range(len(model_i.mol)): 3539 if model_i.mol[j].mol_name != self.structural_data[0].mol[j].mol_name: 3540 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))

Source Code for Module lib.structure.internal.object