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