1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 """Module for the manipulation of the interatomic data structures in the relax data store."""
24
25
26 from copy import deepcopy
27 from numpy import float64, zeros
28 from numpy.linalg import norm
29 from re import search
30 import sys
31 from warnings import warn
32
33
34 from lib.arg_check import is_float
35 from lib.errors import RelaxError, RelaxInteratomInconsistentError, RelaxNoInteratomError, RelaxNoSpinError
36 from lib.io import extract_data, strip, write_data
37 from lib.warnings import RelaxWarning, RelaxZeroVectorWarning
38 from pipe_control import pipes
39 from pipe_control.mol_res_spin import Selection, count_spins, exists_mol_res_spin_data, generate_spin_id_unique, return_spin, spin_loop
40
41
42 -def copy(pipe_from=None, pipe_to=None, spin_id1=None, spin_id2=None, verbose=True):
43 """Copy the interatomic data from one data pipe to another.
44
45 @keyword pipe_from: The data pipe to copy the interatomic data from. This defaults to the current data pipe.
46 @type pipe_from: str
47 @keyword pipe_to: The data pipe to copy the interatomic data to. This defaults to the current data pipe.
48 @type pipe_to: str
49 @keyword spin_id1: The spin ID string of the first atom.
50 @type spin_id1: str
51 @keyword spin_id2: The spin ID string of the second atom.
52 @type spin_id2: str
53 @keyword verbose: A flag which if True will cause info about each spin pair to be printed out.
54 @type verbose: bool
55 """
56
57
58 if pipe_from == None and pipe_to == None:
59 raise RelaxError("The pipe_from and pipe_to arguments cannot both be set to None.")
60 elif pipe_from == None:
61 pipe_from = pipes.cdp_name()
62 elif pipe_to == None:
63 pipe_to = pipes.cdp_name()
64
65
66 pipes.test(pipe_from)
67 pipes.test(pipe_to)
68
69
70 if spin_id1:
71 if count_spins(selection=spin_id1, pipe=pipe_from, skip_desel=False) == 0:
72 raise RelaxNoSpinError(spin_id1, pipe_from)
73 if count_spins(selection=spin_id1, pipe=pipe_to, skip_desel=False) == 0:
74 raise RelaxNoSpinError(spin_id1, pipe_to)
75 if spin_id2:
76 if count_spins(selection=spin_id2, pipe=pipe_from, skip_desel=False) == 0:
77 raise RelaxNoSpinError(spin_id2, pipe_from)
78 if count_spins(selection=spin_id2, pipe=pipe_to, skip_desel=False) == 0:
79 raise RelaxNoSpinError(spin_id2, pipe_to)
80
81
82 if not spin_id1 and not spin_id2:
83 for spin, spin_id in spin_loop(pipe=pipe_from, return_id=True):
84 if not return_spin(spin_id, pipe=pipe_to):
85 raise RelaxNoSpinError(spin_id, pipe_to)
86
87
88 if not exists_data(pipe_from):
89 return
90
91
92 ids = []
93 for interatom in interatomic_loop(selection1=spin_id1, selection2=spin_id2, pipe=pipe_from):
94
95 new_interatom = create_interatom(spin_id1=interatom.spin_id1, spin_id2=interatom.spin_id2, pipe=pipe_to)
96
97
98 for name in dir(interatom):
99
100 if search('^_', name):
101 continue
102
103
104 if name in ['spin_id1', 'spin_id2']:
105 continue
106
107
108 if name in list(interatom.__class__.__dict__.keys()):
109 continue
110
111
112 obj = deepcopy(getattr(interatom, name))
113 setattr(new_interatom, name, obj)
114
115
116 ids.append([repr(interatom.spin_id1), repr(interatom.spin_id2)])
117
118
119 if verbose:
120 write_data(out=sys.stdout, headings=["Spin_ID_1", "Spin_ID_2"], data=ids)
121
122
152
153
155 """Create and return the interatomic data container for the two spins.
156
157 @keyword spin_id1: The spin ID string of the first atom.
158 @type spin_id1: str
159 @keyword spin_id2: The spin ID string of the second atom.
160 @type spin_id2: str
161 @keyword pipe: The data pipe to create the interatomic data container for. This defaults to the current data pipe if not supplied.
162 @type pipe: str or None
163 @keyword verbose: A flag which if True will result printouts.
164 @type verbose: bool
165 @return: The newly created interatomic data container.
166 @rtype: data.interatomic.InteratomContainer instance
167 """
168
169
170 if verbose:
171 print("Creating an interatomic data container between the spins '%s' and '%s'." % (spin_id1, spin_id2))
172
173
174 if pipe == None:
175 pipe = pipes.cdp_name()
176
177
178 dp = pipes.get_pipe(pipe)
179
180
181 spin = return_spin(spin_id1, pipe)
182 if spin == None:
183 raise RelaxNoSpinError(spin_id1)
184 spin = return_spin(spin_id2, pipe)
185 if spin == None:
186 raise RelaxNoSpinError(spin_id2)
187
188
189 for i in range(len(dp.interatomic)):
190 if id_match(spin_id=spin_id1, interatom=dp.interatomic[i], pipe=pipe) and id_match(spin_id=spin_id2, interatom=dp.interatomic[i], pipe=pipe):
191 raise RelaxError("The spin pair %s and %s have already been added." % (spin_id1, spin_id2))
192
193
194 return dp.interatomic.add_item(spin_id1=spin_id1, spin_id2=spin_id2)
195
196
197 -def define(spin_id1=None, spin_id2=None, pipe=None, direct_bond=False, verbose=True):
198 """Set up the magnetic dipole-dipole interaction.
199
200 @keyword spin_id1: The spin identifier string of the first spin of the pair.
201 @type spin_id1: str
202 @keyword spin_id2: The spin identifier string of the second spin of the pair.
203 @type spin_id2: str
204 @param pipe: The data pipe to operate on. Defaults to the current data pipe.
205 @type pipe: str
206 @keyword direct_bond: A flag specifying if the two spins are directly bonded.
207 @type direct_bond: bool
208 @keyword verbose: A flag which if True will result in printouts of the created interatomoic data containers.
209 @type verbose: bool
210 """
211
212
213 if pipe == None:
214 pipe = pipes.cdp_name()
215
216
217 dp = pipes.get_pipe(pipe)
218
219
220 ids = []
221
222
223 if hasattr(dp, 'structure'):
224
225 for mol_name1, res_num1, res_name1, atom_num1, atom_name1, mol_index1, atom_index1 in dp.structure.atom_loop(atom_id=spin_id1, model_num=1, mol_name_flag=True, res_num_flag=True, res_name_flag=True, atom_num_flag=True, atom_name_flag=True, mol_index_flag=True, index_flag=True):
226
227 id1 = generate_spin_id_unique(pipe_cont=dp, mol_name=mol_name1, res_num=res_num1, res_name=res_name1, spin_num=atom_num1, spin_name=atom_name1)
228
229
230 if not return_spin(id1):
231 continue
232
233
234 for mol_name2, res_num2, res_name2, atom_num2, atom_name2, mol_index2, atom_index2 in dp.structure.atom_loop(atom_id=spin_id2, model_num=1, mol_name_flag=True, res_num_flag=True, res_name_flag=True, atom_num_flag=True, atom_name_flag=True, mol_index_flag=True, index_flag=True):
235
236 if direct_bond:
237
238 if mol_name1 != mol_name2:
239 continue
240
241
242 if not dp.structure.are_bonded_index(mol_index1=mol_index1, atom_index1=atom_index1, mol_index2=mol_index2, atom_index2=atom_index2):
243 continue
244
245
246 id2 = generate_spin_id_unique(pipe_cont=dp, mol_name=mol_name2, res_num=res_num2, res_name=res_name2, spin_num=atom_num2, spin_name=atom_name2)
247
248
249 if not return_spin(id2):
250 continue
251
252
253 ids.append([id1, id2])
254
255
256 if ids == []:
257 for spin1, mol_name1, res_num1, res_name1, id1 in spin_loop(spin_id1, pipe=pipe, full_info=True, return_id=True):
258 for spin2, mol_name2, res_num2, res_name2, id2 in spin_loop(spin_id2, pipe=pipe, full_info=True, return_id=True):
259
260 if direct_bond:
261
262 if mol_name1 != mol_name2:
263 continue
264
265
266 if not hasattr(spin1, 'element'):
267 raise RelaxError("The spin '%s' does not have the element type set." % id1)
268 if not hasattr(spin2, 'element'):
269 raise RelaxError("The spin '%s' does not have the element type set." % id2)
270
271
272 pair = False
273 if (spin1.element == 'N' and spin2.element == 'H') or (spin2.element == 'N' and spin1.element == 'H'):
274 pair = True
275 elif (spin1.element == 'C' and spin2.element == 'H') or (spin2.element == 'C' and spin1.element == 'H'):
276 pair = True
277
278
279 if pair and res_num1 != None and res_num1 != res_num2:
280 continue
281 elif pair and res_num1 == None and res_name1 != res_name2:
282 continue
283
284
285 ids.append([id1, id2])
286
287
288 if not len(ids):
289
290 count1 = 0
291 count2 = 0
292 for spin in spin_loop(spin_id1):
293 count1 += 1
294 for spin in spin_loop(spin_id2):
295 count2 += 1
296
297
298 if count1 == 0 and count2 == 0:
299 raise RelaxError("Neither spin IDs '%s' and '%s' match any spins." % (spin_id1, spin_id2))
300 elif count1 == 0:
301 raise RelaxError("The spin ID '%s' matches no spins." % spin_id1)
302 elif count2 == 0:
303 raise RelaxError("The spin ID '%s' matches no spins." % spin_id2)
304 else:
305 raise RelaxError("Unknown error.")
306
307
308 for id1, id2 in ids:
309
310 interatom = return_interatom(id1, id2, pipe=pipe)
311
312
313 if interatom == None:
314 interatom = create_interatom(spin_id1=id1, spin_id2=id2, pipe=pipe)
315
316
317 if interatom.dipole_pair:
318 raise RelaxError("The magnetic dipole-dipole interaction already exists between the spins '%s' and '%s'." % (id1, id2))
319
320
321 interatom.dipole_pair = True
322
323
324 if verbose:
325
326 for i in range(len(ids)):
327 ids[i][0] = repr(ids[i][0])
328 ids[i][1] = repr(ids[i][1])
329
330
331 print("Interatomic interactions are now defined for the following spins:\n")
332 write_data(out=sys.stdout, headings=["Spin_ID_1", "Spin_ID_2"], data=ids)
333
334
336 """Determine if any interatomic data exists.
337
338 @keyword pipe: The data pipe in which the interatomic data will be checked for.
339 @type pipe: str
340 @return: The answer to the question about the existence of data.
341 @rtype: bool
342 """
343
344
345 if pipe == None:
346 pipe = pipes.cdp_name()
347
348
349 pipes.test(pipe)
350
351
352 dp = pipes.get_pipe(pipe)
353
354
355 if dp.interatomic.is_empty():
356 return False
357
358
359 return True
360
361
362 -def id_match(spin_id=None, interatom=None, pipe=None):
363 """Test if the spin ID matches one of the two spins of the given container.
364
365 @keyword spin_id: The spin ID string of the first atom.
366 @type spin_id: str
367 @keyword interatom: The interatomic data container.
368 @type interatom: InteratomContainer instance
369 @keyword pipe: The data pipe containing the interatomic data container. Defaults to the current data pipe.
370 @type pipe: str or None
371 @return: True if the spin ID matches one of the two spins, False otherwise.
372 @rtype: bool
373 """
374
375
376 spin1 = return_spin(interatom.spin_id1, pipe=pipe)
377 spin2 = return_spin(interatom.spin_id2, pipe=pipe)
378
379
380 if spin1 == None or spin2 == None:
381 return False
382
383
384 if spin_id in spin1._spin_ids or spin_id in spin2._spin_ids:
385 return True
386
387
388 return False
389
390
391 -def interatomic_loop(selection1=None, selection2=None, pipe=None, skip_desel=True):
392 """Generator function for looping over all the interatomic data containers.
393
394 @keyword selection1: The optional spin ID selection of the first atom.
395 @type selection1: str
396 @keyword selection2: The optional spin ID selection of the second atom.
397 @type selection2: str
398 @keyword pipe: The data pipe containing the spin. Defaults to the current data pipe.
399 @type pipe: str
400 @keyword skip_desel: A flag which if True will cause only selected interatomic data containers to be returned.
401 @type skip_desel: bool
402 """
403
404
405 if pipe == None:
406 pipe = pipes.cdp_name()
407
408
409 dp = pipes.get_pipe(pipe)
410
411
412 select_obj = None
413 select_obj1 = None
414 select_obj2 = None
415 if selection1 and selection2:
416 select_obj1 = Selection(selection1)
417 select_obj2 = Selection(selection2)
418 elif selection1:
419 select_obj = Selection(selection1)
420 elif selection2:
421 select_obj = Selection(selection2)
422
423
424 for i in range(len(dp.interatomic)):
425
426 if skip_desel and not dp.interatomic[i].select:
427 continue
428
429
430 interatom = dp.interatomic[i]
431 mol_index1, res_index1, spin_index1 = cdp.mol._spin_id_lookup[interatom.spin_id1]
432 mol_index2, res_index2, spin_index2 = cdp.mol._spin_id_lookup[interatom.spin_id2]
433 mol1 = cdp.mol[mol_index1]
434 res1 = cdp.mol[mol_index1].res[res_index1]
435 spin1 = cdp.mol[mol_index1].res[res_index1].spin[spin_index1]
436 mol2 = cdp.mol[mol_index2]
437 res2 = cdp.mol[mol_index2].res[res_index2]
438 spin2 = cdp.mol[mol_index2].res[res_index2].spin[spin_index2]
439
440
441 if select_obj:
442 sel1 = select_obj.contains_spin(spin_name=spin1.name, spin_num=spin1.num, res_name=res1.name, res_num=res1.num, mol=mol1.name)
443 sel2 = select_obj.contains_spin(spin_name=spin2.name, spin_num=spin2.num, res_name=res2.name, res_num=res2.num, mol=mol2.name)
444 if select_obj1:
445 sel11 = select_obj1.contains_spin(spin_name=spin1.name, spin_num=spin1.num, res_name=res1.name, res_num=res1.num, mol=mol1.name)
446 sel12 = select_obj1.contains_spin(spin_name=spin2.name, spin_num=spin2.num, res_name=res2.name, res_num=res2.num, mol=mol2.name)
447 if select_obj2:
448 sel21 = select_obj2.contains_spin(spin_name=spin1.name, spin_num=spin1.num, res_name=res1.name, res_num=res1.num, mol=mol1.name)
449 sel22 = select_obj2.contains_spin(spin_name=spin2.name, spin_num=spin2.num, res_name=res2.name, res_num=res2.num, mol=mol2.name)
450
451
452 if select_obj:
453 if not sel1 and not sel2:
454 continue
455 if select_obj1:
456 if not (sel11 or sel12):
457 continue
458 if select_obj2:
459 if not (sel21 or sel22):
460 continue
461
462
463 yield interatom
464
465
466 -def read_dist(file=None, dir=None, unit='meter', spin_id1_col=None, spin_id2_col=None, data_col=None, sep=None):
467 """Set up the magnetic dipole-dipole interaction.
468
469 @keyword file: The name of the file to open.
470 @type file: str
471 @keyword dir: The directory containing the file (defaults to the current directory if None).
472 @type dir: str or None
473 @keyword unit: The measurement unit. This can be either 'meter' or 'Angstrom'.
474 @type unit: str
475 @keyword spin_id1_col: The column containing the spin ID strings of the first spin.
476 @type spin_id1_col: int
477 @keyword spin_id2_col: The column containing the spin ID strings of the second spin.
478 @type spin_id2_col: int
479 @keyword data_col: The column containing the averaged distances in meters.
480 @type data_col: int or None
481 @keyword sep: The column separator which, if None, defaults to whitespace.
482 @type sep: str or None
483 """
484
485
486 if unit not in ['meter', 'Angstrom']:
487 raise RelaxError("The measurement unit of '%s' must be one of 'meter' or 'Angstrom'." % unit)
488
489
490 pipes.test()
491
492
493 if not exists_mol_res_spin_data():
494 raise RelaxNoSequenceError
495
496
497 file_data = extract_data(file, dir, sep=sep)
498 file_data = strip(file_data, comments=True)
499
500
501 data = []
502 for line in file_data:
503
504 if spin_id1_col > len(line):
505 warn(RelaxWarning("The data %s is invalid, no first spin ID column can be found." % line))
506 continue
507 if spin_id2_col > len(line):
508 warn(RelaxWarning("The data %s is invalid, no second spin ID column can be found." % line))
509 continue
510 if data_col and data_col > len(line):
511 warn(RelaxWarning("The data %s is invalid, no data column can be found." % line))
512 continue
513
514
515 spin_id1 = line[spin_id1_col-1]
516 spin_id2 = line[spin_id2_col-1]
517 ave_dist = None
518 if data_col:
519 ave_dist = line[data_col-1]
520
521
522 if ave_dist != None:
523 try:
524 ave_dist = float(ave_dist)
525 except ValueError:
526 warn(RelaxWarning("The averaged distance of '%s' from the line %s is invalid." % (ave_dist, line)))
527 continue
528
529
530 if unit == 'Angstrom':
531 ave_dist = ave_dist * 1e-10
532
533
534 interatom = return_interatom(spin_id1, spin_id2)
535
536
537 if interatom == None:
538 interatom = create_interatom(spin_id1=spin_id1, spin_id2=spin_id2, verbose=True)
539
540
541 interatom.r = ave_dist
542
543
544 data.append([repr(interatom.spin_id1), repr(interatom.spin_id2), repr(ave_dist)])
545
546
547 if not len(data):
548 raise RelaxError("No data could be extracted from the file.")
549
550
551 print("The following averaged distances have been read:\n")
552 write_data(out=sys.stdout, headings=["Spin_ID_1", "Spin_ID_2", "Ave_distance(meters)"], data=data)
553
554
556 """Return the list of interatomic data containers for the two spins.
557
558 @keyword spin_id1: The spin ID string of the first atom.
559 @type spin_id1: str
560 @keyword spin_id2: The spin ID string of the second atom.
561 @type spin_id2: str
562 @keyword pipe: The data pipe holding the container. Defaults to the current data pipe.
563 @type pipe: str or None
564 @return: The matching interatomic data container, if it exists.
565 @rtype: data.interatomic.InteratomContainer instance or None
566 """
567
568
569 if spin_id1 == None:
570 raise RelaxError("The first spin ID must be supplied.")
571 if spin_id2 == None:
572 raise RelaxError("The second spin ID must be supplied.")
573
574
575 if pipe == None:
576 pipe = pipes.cdp_name()
577
578
579 dp = pipes.get_pipe(pipe)
580
581
582 for i in range(len(dp.interatomic)):
583 if id_match(spin_id=spin_id1, interatom=dp.interatomic[i], pipe=pipe) and id_match(spin_id=spin_id2, interatom=dp.interatomic[i], pipe=pipe):
584 return dp.interatomic[i]
585
586
587 return None
588
589
622
623
624 -def set_dist(spin_id1=None, spin_id2=None, ave_dist=None, unit='meter'):
625 """Set up the magnetic dipole-dipole interaction.
626
627 @keyword spin_id1: The spin identifier string of the first spin of the pair.
628 @type spin_id1: str
629 @keyword spin_id2: The spin identifier string of the second spin of the pair.
630 @type spin_id2: str
631 @keyword ave_dist: The r^-3 averaged interatomic distance.
632 @type ave_dist: float
633 @keyword unit: The measurement unit. This can be either 'meter' or 'Angstrom'.
634 @type unit: str
635 """
636
637
638 if unit not in ['meter', 'Angstrom']:
639 raise RelaxError("The measurement unit of '%s' must be one of 'meter' or 'Angstrom'." % unit)
640
641
642 if unit == 'Angstrom':
643 ave_dist = ave_dist * 1e-10
644
645
646 sel_obj1 = Selection(spin_id1)
647 sel_obj2 = Selection(spin_id2)
648
649
650 data = []
651 for interatom in interatomic_loop():
652
653 mol_name1, res_num1, res_name1, spin1 = return_spin(interatom.spin_id1, full_info=True)
654 mol_name2, res_num2, res_name2, spin2 = return_spin(interatom.spin_id2, full_info=True)
655
656
657 if not (sel_obj1.contains_spin(spin_num=spin1.num, spin_name=spin1.name, res_num=res_num1, res_name=res_name1, mol=mol_name1) and sel_obj2.contains_spin(spin_num=spin2.num, spin_name=spin2.name, res_num=res_num2, res_name=res_name2, mol=mol_name2)) and not (sel_obj2.contains_spin(spin_num=spin1.num, spin_name=spin1.name, res_num=res_num1, res_name=res_name1, mol=mol_name1) and sel_obj1.contains_spin(spin_num=spin2.num, spin_name=spin2.name, res_num=res_num2, res_name=res_name2, mol=mol_name2)):
658 continue
659
660
661 interatom.r = ave_dist
662
663
664 data.append([repr(interatom.spin_id1), repr(interatom.spin_id2), repr(ave_dist)])
665
666
667 if not len(data):
668 raise RelaxError("No data could be set.")
669
670
671 print("The following averaged distances have been set:\n")
672 write_data(out=sys.stdout, headings=["Spin_ID_1", "Spin_ID_2", "Ave_distance(meters)"], data=data)
673
674
676 """Extract the bond vectors from the loaded structures and store them in the spin container.
677
678 @keyword ave: A flag which if True will cause the average of all vectors to be calculated.
679 @type ave: bool
680 """
681
682
683 pipes.test()
684
685
686 if not exists_data():
687 raise RelaxNoInteratomError
688
689
690 if ave:
691 print("Averaging all vectors.")
692 else:
693 print("No averaging of the vectors.")
694
695
696 no_vectors = True
697 pos_info = False
698 for interatom in interatomic_loop(skip_desel=False):
699
700 spin1 = return_spin(interatom.spin_id1)
701 spin2 = return_spin(interatom.spin_id2)
702
703
704 if not hasattr(spin1, 'pos'):
705 continue
706 if not hasattr(spin2, 'pos'):
707 continue
708
709
710 pos_info = True
711
712
713 if is_float(spin1.pos[0], raise_error=False) and is_float(spin2.pos[0], raise_error=False):
714
715 vector_list = [spin2.pos - spin1.pos]
716
717
718 elif is_float(spin1.pos[0], raise_error=False) or is_float(spin2.pos[0], raise_error=False):
719
720 if is_float(spin2.pos[0], raise_error=False):
721 vector_list = []
722 for i in range(len(spin1.pos)):
723 vector_list.append(spin2.pos - spin1.pos[i])
724
725
726 else:
727 vector_list = []
728 for i in range(len(spin2.pos)):
729 vector_list.append(spin2.pos[i] - spin1.pos)
730
731
732 else:
733
734 if len(spin1.pos) != len(spin2.pos):
735 raise RelaxError("The spin '%s' consists of %s positions whereas the spin '%s' consists of %s - these numbers must match." % (interatom.spin_id1, len(spin1.pos), interatom.spin_id1, len(spin1.pos)))
736
737
738 vector_list = []
739 for i in range(len(spin1.pos)):
740 vector_list.append(spin2.pos[i] - spin1.pos[i])
741
742
743 for i in range(len(vector_list)):
744
745 norm_factor = norm(vector_list[i])
746
747
748 if norm_factor == 0.0:
749 warn(RelaxZeroVectorWarning(spin_id1=interatom.spin_id1, spin_id2=interatom.spin_id2))
750
751
752 else:
753 vector_list[i] = vector_list[i] / norm_factor
754
755
756 if ave:
757 ave_vector = zeros(3, float64)
758 for i in range(len(vector_list)):
759 ave_vector = ave_vector + vector_list[i]
760 vector_list = [ave_vector / len(vector_list)]
761
762
763 if len(vector_list) == 1:
764 vector_list = vector_list[0]
765
766
767 setattr(interatom, 'vector', vector_list)
768
769
770 no_vectors = False
771
772
773 num = 1
774 if not is_float(vector_list[0], raise_error=False):
775 num = len(vector_list)
776 plural = 's'
777 if num == 1:
778 plural = ''
779 if spin1.name:
780 spin1_str = spin1.name
781 else:
782 spin1_str = spin1.num
783 if spin2.name:
784 spin2_str = spin2.name
785 else:
786 spin2_str = spin2.num
787 print("Calculated %s %s-%s unit vector%s between the spins '%s' and '%s'." % (num, spin1_str, spin2_str, plural, interatom.spin_id1, interatom.spin_id2))
788
789
790 if not pos_info:
791 raise RelaxError("Positional information could not be found for any spins.")
792
793
794 if no_vectors:
795 raise RelaxError("No vectors could be extracted.")
796