1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24 """Module for the manipulation of pseudo-contact shift data."""
25
26
27 from copy import deepcopy
28 from math import ceil, floor, pi, sqrt
29 from numpy import array, eye, float64, int32, ones, std, zeros
30 from numpy.linalg import norm
31 from numpy.random import multivariate_normal
32 import sys
33 from warnings import warn
34
35
36 from lib.alignment.pcs import ave_pcs_tensor, pcs_tensor
37 from lib.check_types import is_float
38 from lib.errors import RelaxError, RelaxNoAlignError, RelaxNoPdbError, RelaxNoPCSError, RelaxNoSequenceError
39 from lib.io import open_write_file, write_data
40 from lib.periodic_table import periodic_table
41 from lib.physical_constants import pcs_constant
42 from lib.plotting.api import write_xy_data, write_xy_header
43 from lib.sequence import read_spin_data, write_spin_data
44 from lib.warnings import RelaxWarning, RelaxNoSpinWarning
45 from pipe_control import pipes
46 from pipe_control.align_tensor import get_tensor_index, get_tensor_object, opt_uses_align_data, opt_uses_tensor
47 from pipe_control.mol_res_spin import exists_mol_res_spin_data, generate_spin_id_unique, is_pseudoatom, return_spin, spin_loop
48 from pipe_control.pipes import check_pipe
49
50
52 """Back calculate the PCS from the alignment tensor.
53
54 @keyword align_id: The alignment tensor ID string.
55 @type align_id: str
56 """
57
58
59 check_pipe_setup(pcs_id=align_id, sequence=True, N=True, tensors=True, paramag_centre=True)
60
61
62 if align_id:
63 align_ids = [align_id]
64 else:
65 align_ids = cdp.align_ids
66
67
68 for align_id in align_ids:
69
70 if not hasattr(cdp, 'pcs_ids'):
71 cdp.pcs_ids = []
72
73
74 if align_id not in cdp.pcs_ids:
75 cdp.pcs_ids.append(align_id)
76
77
78 weights = ones(cdp.N, float64) / cdp.N
79
80
81 unit_vect = zeros((cdp.N, 3), float64)
82
83
84 count = 0
85 for spin in spin_loop():
86
87 if not hasattr(spin, 'pos'):
88 continue
89
90
91 pos = spin.pos
92 if type(pos[0]) in [float, float64]:
93 pos = [pos] * cdp.N
94
95
96 for id in align_ids:
97
98 vect = zeros((cdp.N, 3), float64)
99 r = zeros(cdp.N, float64)
100 dj = zeros(cdp.N, float64)
101 for c in range(cdp.N):
102
103 vect[c] = pos[c] - cdp.paramagnetic_centre
104
105
106 r[c] = norm(vect[c])
107
108
109 if r[c]:
110 vect[c] = vect[c] / r[c]
111
112
113 dj[c] = pcs_constant(cdp.temperature[id], cdp.spectrometer_frq[id] * 2.0 * pi / periodic_table.gyromagnetic_ratio('1H'), r[c]/1e10)
114
115
116 if not hasattr(spin, 'pcs_bc'):
117 spin.pcs_bc = {}
118
119
120 spin.pcs_bc[id] = ave_pcs_tensor(dj, vect, cdp.N, cdp.align_tensors[get_tensor_index(align_id=id)].A, weights=weights) * 1e6
121
122
123 count += 1
124
125
126 if not count:
127 warn(RelaxWarning("No PCSs have been back calculated, probably due to missing spin position information."))
128
129
130 -def centre(pos=None, atom_id=None, pipe=None, verbosity=1, ave_pos=False, force=False):
131 """Specify the atom in the loaded structure corresponding to the paramagnetic centre.
132
133 @keyword pos: The atomic position. If set, the atom_id string will be ignored.
134 @type pos: list of float
135 @keyword atom_id: The atom identification string.
136 @type atom_id: str
137 @keyword pipe: An alternative data pipe to extract the paramagnetic centre from.
138 @type pipe: None or str
139 @keyword verbosity: The amount of information to print out. The bigger the number, the more information.
140 @type verbosity: int
141 @keyword ave_pos: A flag which if True causes the atomic positions from multiple models to be averaged.
142 @type ave_pos: bool
143 @keyword force: A flag which if True will cause the current PCS centre to be overwritten.
144 """
145
146
147 if pipe == None:
148 pipe = pipes.cdp_name()
149
150
151 check_pipe_setup(pipe=pipe)
152
153
154 source_dp = pipes.get_pipe(pipe)
155
156
157 if not hasattr(source_dp, 'structure'):
158 raise RelaxNoPdbError
159
160
161 if not force and hasattr(cdp, 'paramagnetic_centre'):
162 raise RelaxError("The paramagnetic centre has already been set to the coordinates " + repr(cdp.paramagnetic_centre) + ".")
163
164
165 if pos != None:
166 centre = array(pos)
167 num_pos = 1
168 full_pos_list = []
169
170
171 else:
172
173 centre = zeros(3, float64)
174 full_pos_list = []
175 num_pos = 0
176 for spin, spin_id in spin_loop(atom_id, pipe=pipe, return_id=True):
177
178 if not hasattr(spin, 'pos'):
179 continue
180
181
182 if isinstance(spin.pos[0], float) or isinstance(spin.pos[0], float64):
183 pos_list = [spin.pos]
184 else:
185 pos_list = spin.pos
186
187
188 for pos in pos_list:
189 full_pos_list.append(pos)
190 centre = centre + array(pos)
191 num_pos = num_pos + 1
192
193
194 if not num_pos:
195 raise RelaxError("No positional information could be found for the spin '%s'." % atom_id)
196
197
198 centre = centre / float(num_pos)
199
200
201 if verbosity:
202 print("Paramagnetic centres located at:")
203 for pos in full_pos_list:
204 print(" [%8.3f, %8.3f, %8.3f]" % (pos[0], pos[1], pos[2]))
205 print("\nAverage paramagnetic centre located at:")
206 print(" [%8.3f, %8.3f, %8.3f]" % (centre[0], centre[1], centre[2]))
207
208
209 if ave_pos:
210 if verbosity:
211 print("\nUsing the average paramagnetic position.")
212 cdp.paramagnetic_centre = centre
213 else:
214 if verbosity:
215 print("\nUsing all paramagnetic positions.")
216 cdp.paramagnetic_centre = full_pos_list
217
218
219 -def check_pipe_setup(pipe=None, pcs_id=None, sequence=False, N=False, tensors=False, pcs=False, paramag_centre=False):
220 """Check that the current data pipe has been setup sufficiently.
221
222 @keyword pipe: The data pipe to check, defaulting to the current pipe.
223 @type pipe: None or str
224 @keyword pcs_id: The PCS ID string to check for in cdp.pcs_ids.
225 @type pcs_id: None or str
226 @keyword sequence: A flag which when True will invoke the sequence data check.
227 @type sequence: bool
228 @keyword N: A flag which if True will check that cdp.N is set.
229 @type N: bool
230 @keyword tensors: A flag which if True will check that alignment tensors exist.
231 @type tensors: bool
232 @keyword pcs: A flag which if True will check that PCSs exist.
233 @type pcs: bool
234 @keyword paramag_centre: A flag which if True will check that the paramagnetic centre has been set.
235 @type paramag_centre: bool
236 """
237
238
239 if pipe == None:
240 pipe = pipes.cdp_name()
241
242
243 dp = pipes.get_pipe(pipe)
244
245
246 check_pipe(pipe)
247
248
249 if sequence and not exists_mol_res_spin_data(pipe):
250 raise RelaxNoSequenceError(pipe)
251
252
253 if N and not hasattr(dp, 'N'):
254 raise RelaxError("The number of states N has not been set.")
255
256
257 if tensors and (not hasattr(dp, 'align_tensors') or len(dp.align_tensors) == 0):
258 raise RelaxNoTensorError('alignment')
259
260
261 if pcs_id and (not hasattr(dp, 'align_ids') or pcs_id not in dp.align_ids):
262 raise RelaxNoAlignError(pcs_id, pipe)
263
264
265 if pcs and not hasattr(dp, 'align_ids'):
266 raise RelaxNoAlignError()
267 if pcs and not hasattr(dp, 'pcs_ids'):
268 raise RelaxNoPCSError()
269 elif pcs and pcs_id and pcs_id not in dp.pcs_ids:
270 raise RelaxNoPCSError(pcs_id)
271
272
273 if paramag_centre and not hasattr(cdp, 'paramagnetic_centre'):
274 raise RelaxError("The paramagnetic centre has not been defined.")
275
276
277 -def copy(pipe_from=None, pipe_to=None, align_id=None, back_calc=True):
278 """Copy the PCS data from one data pipe to another.
279
280 @keyword pipe_from: The data pipe to copy the PCS data from. This defaults to the current data pipe.
281 @type pipe_from: str
282 @keyword pipe_to: The data pipe to copy the PCS data to. This defaults to the current data pipe.
283 @type pipe_to: str
284 @keyword align_id: The alignment ID string.
285 @type align_id: str
286 @keyword back_calc: A flag which if True will cause any back-calculated RDCs present to also be copied with the real values and errors.
287 @type back_calc: bool
288 """
289
290
291 if pipe_from == None and pipe_to == None:
292 raise RelaxError("The pipe_from and pipe_to arguments cannot both be set to None.")
293 elif pipe_from == None:
294 pipe_from = pipes.cdp_name()
295 elif pipe_to == None:
296 pipe_to = pipes.cdp_name()
297
298
299 check_pipe_setup(pipe=pipe_from, pcs_id=align_id, sequence=True, pcs=True)
300 check_pipe_setup(pipe=pipe_to, sequence=True)
301
302
303 dp_from = pipes.get_pipe(pipe_from)
304 dp_to = pipes.get_pipe(pipe_to)
305
306
307 if align_id == None:
308 align_ids = dp_from.align_ids
309 else:
310 align_ids = [align_id]
311
312
313 if not hasattr(dp_to, 'align_ids'):
314 dp_to.align_ids = []
315 if not hasattr(dp_to, 'pcs_ids'):
316 dp_to.pcs_ids = []
317
318
319 for align_id in align_ids:
320
321 print("\nCoping PCSs for the alignment ID '%s'." % align_id)
322
323
324 if align_id not in dp_to.align_ids and align_id not in dp_to.align_ids:
325 dp_to.align_ids.append(align_id)
326 if align_id in dp_from.pcs_ids and align_id not in dp_to.pcs_ids:
327 dp_to.pcs_ids.append(align_id)
328
329
330 data = []
331 for spin_from, spin_id in spin_loop(return_id=True, skip_desel=True, pipe=pipe_from):
332
333 spin_to = return_spin(spin_id=spin_id, pipe=pipe_to)
334
335
336 if spin_to == None:
337 warn(RelaxWarning("The spin container for the spin '%s' cannot be found in the target data pipe." % spin_id))
338 continue
339
340
341 if (not hasattr(spin_from, 'pcs') or not align_id in spin_from.pcs) and (not hasattr(spin_from, 'pcs_err') or not align_id in spin_from.pcs_err):
342 continue
343
344
345 if hasattr(spin_from, 'pcs') and not hasattr(spin_to, 'pcs'):
346 spin_to.pcs = {}
347 if back_calc and hasattr(spin_from, 'pcs_bc') and not hasattr(spin_to, 'pcs_bc'):
348 spin_to.pcs_bc = {}
349 if hasattr(spin_from, 'pcs_err') and not hasattr(spin_to, 'pcs_err'):
350 spin_to.pcs_err = {}
351
352
353 value = None
354 error = None
355 value_bc = None
356 if hasattr(spin_from, 'pcs'):
357 value = spin_from.pcs[align_id]
358 spin_to.pcs[align_id] = value
359 if back_calc and hasattr(spin_from, 'pcs_bc'):
360 value_bc = spin_from.pcs_bc[align_id]
361 spin_to.pcs_bc[align_id] = value_bc
362 if hasattr(spin_from, 'pcs_err'):
363 error = spin_from.pcs_err[align_id]
364 spin_to.pcs_err[align_id] = error
365
366
367 data.append([spin_id])
368 if is_float(value):
369 data[-1].append("%20.15f" % value)
370 else:
371 data[-1].append("%20s" % value)
372 if back_calc:
373 if is_float(value_bc):
374 data[-1].append("%20.15f" % value_bc)
375 else:
376 data[-1].append("%20s" % value_bc)
377 if is_float(error):
378 data[-1].append("%20.15f" % error)
379 else:
380 data[-1].append("%20s" % error)
381
382
383 print("The following PCSs have been copied:\n")
384 if back_calc:
385 write_data(out=sys.stdout, headings=["Spin_ID", "Value", "Back-calculated", "Error"], data=data)
386 else:
387 write_data(out=sys.stdout, headings=["Spin_ID", "Value", "Error"], data=data)
388
389
390 -def corr_plot(format=None, title=None, subtitle=None, file=None, dir=None, force=False):
391 """Generate a correlation plot of the measured vs. back-calculated PCSs.
392
393 @keyword format: The format for the plot file. The following values are accepted: 'grace', a Grace plot; None, a plain text file.
394 @type format: str or None
395 @keyword title: The title for the plot, overriding the default.
396 @type title: None or str
397 @keyword subtitle: The subtitle for the plot, overriding the default.
398 @type subtitle: None or str
399 @keyword file: The file name or object to write to.
400 @type file: str or file object
401 @keyword dir: The name of the directory to place the file into (defaults to the current directory).
402 @type dir: str
403 @keyword force: A flag which if True will cause any pre-existing file to be overwritten.
404 @type force: bool
405 """
406
407
408 check_pipe_setup(sequence=True)
409
410
411 if not hasattr(cdp, 'pcs_ids') or not cdp.pcs_ids:
412 warn(RelaxWarning("No PCS data exists, skipping file creation."))
413 return
414
415
416 file = open_write_file(file, dir, force)
417
418
419 data = []
420 orig_title = title
421 if orig_title == None:
422 title = "PCS correlation plot"
423 axis_labels = ["Back-calculated PCS (ppm)", "Measured PCS (ppm)"]
424
425
426 data.append([[-100, -100, 0], [100, 100, 0]])
427
428
429 types = []
430 for spin in spin_loop():
431 if not hasattr(spin, 'element'):
432 if None not in types:
433 types.append(None)
434 elif spin.element not in types:
435 types.append(spin.element)
436
437
438 min_pcs = 1e100
439 max_pcs = -1e100
440 for align_id in cdp.pcs_ids:
441
442 for i in range(len(types)):
443
444 data.append([])
445
446
447 err_flag = False
448 for spin in spin_loop():
449
450 if not spin.select:
451 continue
452
453
454 if hasattr(spin, 'pcs_err') and align_id in spin.pcs_err:
455 err_flag = True
456 break
457
458
459 for spin, spin_id in spin_loop(return_id=True):
460
461 if not spin.select:
462 continue
463
464
465 if hasattr(spin, 'element') and spin.element != types[i]:
466 continue
467
468
469 if not hasattr(spin, 'pcs') or not hasattr(spin, 'pcs_bc') or not align_id in spin.pcs or not align_id in spin.pcs_bc:
470 continue
471 if spin.pcs[align_id] == None or spin.pcs_bc[align_id] == None:
472 continue
473
474
475 data[-1].append([spin.pcs_bc[align_id], spin.pcs[align_id]])
476
477
478 if spin.pcs[align_id] < min_pcs:
479 min_pcs = spin.pcs[align_id]
480 if spin.pcs_bc[align_id] < min_pcs:
481 min_pcs = spin.pcs_bc[align_id]
482
483
484 if spin.pcs[align_id] > max_pcs:
485 max_pcs = spin.pcs[align_id]
486 if spin.pcs_bc[align_id] > max_pcs:
487 max_pcs = spin.pcs_bc[align_id]
488
489
490 if err_flag:
491 if hasattr(spin, 'pcs_err') and align_id in spin.pcs_err:
492 data[-1][-1].append(spin.pcs_err[align_id])
493 else:
494 data[-1][-1].append(None)
495
496
497 data[-1][-1].append(spin_id)
498
499
500 size = len(data)
501
502
503 max_pcs = float(ceil(max_pcs))
504 min_pcs = float(floor(min_pcs))
505
506
507 data = [data]
508
509
510 if err_flag:
511 graph_type = 'xydy'
512 else:
513 graph_type = 'xy'
514
515
516 if format == 'grace':
517
518 set_names = [None]
519 for i in range(len(cdp.pcs_ids)):
520 for j in range(len(types)):
521 set_names.append("%s (%s)" % (cdp.pcs_ids[i], types[j]))
522
523
524 write_xy_header(format=format, file=file, title=title, subtitle=subtitle, world=[[min_pcs, min_pcs, max_pcs, max_pcs]], sets=[size], set_names=[set_names], linestyle=[[2]+[0]*size], data_type=['pcs_bc', 'pcs'], axis_labels=[axis_labels], tick_major_spacing=[[1, 1]], tick_minor_count=[[9, 9]], legend_pos=[[1, 0.5]])
525
526
527 write_xy_data(format=format, data=data, file=file, graph_type=graph_type, autoscale=False)
528
529
531 """Delete the PCS data corresponding to the alignment ID.
532
533 @keyword align_id: The alignment tensor ID string. If not specified, all data will be deleted.
534 @type align_id: str or None
535 """
536
537
538 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True)
539
540
541 if not align_id:
542 ids = deepcopy(cdp.pcs_ids)
543 else:
544 ids = [align_id]
545
546
547 for id in ids:
548
549 cdp.pcs_ids.pop(cdp.pcs_ids.index(id))
550
551
552 if hasattr(cdp, 'pcs_data_types') and id in cdp.pcs_data_types:
553 cdp.pcs_data_types.pop(id)
554
555
556 for spin in spin_loop():
557
558 if hasattr(spin, 'pcs') and id in spin.pcs:
559 spin.pcs.pop(id)
560
561
562 if hasattr(spin, 'pcs_err') and id in spin.pcs_err:
563 spin.pcs_err.pop(id)
564
565
566 if not hasattr(cdp, 'rdc_ids') or id not in cdp.rdc_ids:
567 cdp.align_ids.pop(cdp.align_ids.index(id))
568
569
570 -def display(align_id=None, bc=False):
571 """Display the PCS data corresponding to the alignment ID.
572
573 @keyword align_id: The alignment tensor ID string.
574 @type align_id: str
575 @keyword bc: The back-calculation flag which if True will cause the back-calculated rather than measured data to be displayed.
576 @type bc: bool
577 """
578
579
580 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True)
581
582
583 write(align_id=align_id, file=sys.stdout, bc=bc)
584
585
587 """Determine if the PCS data for the given alignment ID is needed for optimisation.
588
589 @param align_id: The alignment ID string.
590 @type align_id: str
591 @return: True if the PCS data is to be used for optimisation, False otherwise.
592 @rtype: bool
593 """
594
595
596 if not hasattr(cdp, 'pcs_ids'):
597 return False
598
599
600 if align_id not in cdp.pcs_ids:
601 return False
602
603
604 tensor_flag = opt_uses_tensor(get_tensor_object(align_id))
605
606
607 pos_flag = False
608 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed:
609 pos_flag = True
610
611
612 prob_flag = False
613 if cdp.model == 'population':
614 prob_flag = True
615
616
617 if not tensor_flag and not pos_flag and not prob_flag:
618 return False
619
620
621 return True
622
623
624 -def q_factors(spin_id=None, sim_index=None, verbosity=1):
625 """Calculate the Q factors for the PCS data.
626
627 @keyword spin_id: The spin ID string used to restrict the Q factor calculation to a subset of all spins.
628 @type spin_id: None or str
629 @keyword sim_index: The optional Monte Carlo simulation index.
630 @type sim_index: None or int
631 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity.
632 @type verbosity: int
633 """
634
635
636 check_pipe_setup(sequence=True)
637
638
639 sim_flag = (sim_index != None)
640
641
642 if not hasattr(cdp, 'pcs_ids') or not len(cdp.pcs_ids):
643 warn(RelaxWarning("No PCS data exists, Q factors cannot be calculated."))
644 return
645
646
647 if not sim_flag:
648 cdp.q_factors_pcs_norm_squared_sum = {}
649 else:
650 if not hasattr(cdp, 'q_factors_pcs_norm_squared_sum_sim'):
651 cdp.q_factors_pcs_norm_squared_sum_sim = {}
652
653
654 for align_id in cdp.pcs_ids:
655
656 pcs2_sum = 0.0
657 sse = 0.0
658 if sim_flag and align_id not in cdp.q_factors_pcs_norm_squared_sum_sim:
659 cdp.q_factors_pcs_norm_squared_sum_sim[align_id] = [None] * cdp.sim_number
660
661
662 spin_count = 0
663 pcs_data = False
664 pcs_bc_data = False
665 for spin in spin_loop(spin_id):
666 pcs_data_align_id = False
667 pcs_bc_data_align_id = False
668
669
670 if not spin.select:
671 continue
672
673
674 spin_count += 1
675
676
677 if not sim_flag:
678 if hasattr(spin, 'pcs') and align_id in spin.pcs and spin.pcs[align_id] != None:
679 pcs_data_align_id = True
680 if hasattr(spin, 'pcs_bc') and align_id in spin.pcs_bc and spin.pcs_bc[align_id] != None:
681 pcs_bc_data_align_id = True
682 else:
683 if hasattr(spin, 'pcs_sim') and align_id in spin.pcs_sim and spin.pcs_sim[align_id][sim_index] != None:
684 pcs_data_align_id = True
685 if hasattr(spin, 'pcs_sim_bc') and align_id in spin.pcs_sim_bc and spin.pcs_sim_bc[align_id][sim_index] != None:
686 pcs_bc_data_align_id = True
687
688
689 pcs_data = pcs_data or pcs_data_align_id
690 pcs_bc_data = pcs_bc_data or pcs_bc_data_align_id
691
692
693 if not pcs_data_align_id or not pcs_bc_data_align_id:
694 continue
695
696
697 if not sim_flag:
698 pcs = spin.pcs[align_id]
699 pcs_bc = spin.pcs_bc[align_id]
700 else:
701 pcs = spin.pcs_sim[align_id][sim_index]
702 pcs_bc = spin.pcs_sim_bc[align_id][sim_index]
703
704
705 sse = sse + (pcs - pcs_bc)**2
706
707
708 pcs2_sum = pcs2_sum + pcs**2
709
710
711 if pcs2_sum:
712 Q = sqrt(sse / pcs2_sum)
713 if sim_flag:
714 cdp.q_factors_pcs_norm_squared_sum_sim[align_id][sim_index] = Q
715 else:
716 cdp.q_factors_pcs_norm_squared_sum[align_id] = Q
717
718
719 if not spin_count:
720 if not sim_flag:
721 warn(RelaxWarning("No spins have been used in the calculation, skipping the PCS Q factor calculation."))
722 return
723 if not pcs_data:
724 if not sim_flag:
725 warn(RelaxWarning("No PCS data can be found for the alignment ID '%s', skipping the PCS Q factor calculation for this alignment." % align_id))
726 continue
727 if not pcs_bc_data:
728 if not sim_flag:
729 warn(RelaxWarning("No back-calculated PCS data can be found for the alignment ID '%s', skipping the PCS Q factor calculation for this alignment." % align_id))
730 continue
731
732
733 if verbosity:
734 print("\nPCS Q factors normalised by the sum of PCSs squared:")
735 for align_id in cdp.pcs_ids:
736 if align_id in cdp.q_factors_pcs_norm_squared_sum:
737 if sim_flag:
738 print(" Alignment ID '%s': %.3f" % (align_id, cdp.q_factors_pcs_norm_squared_sum_sim[align_id][sim_index]))
739 else:
740 print(" Alignment ID '%s': %.3f" % (align_id, cdp.q_factors_pcs_norm_squared_sum[align_id]))
741
742
743 if sim_flag:
744 if not hasattr(cdp, 'q_pcs_norm_squared_sum_sim'):
745 cdp.q_pcs_norm_squared_sum_sim = [None] * cdp.sim_number
746 cdp.q_pcs_norm_squared_sum_sim[sim_index] = 0.0
747 for id in cdp.q_factors_pcs_norm_squared_sum_sim:
748 cdp.q_pcs_norm_squared_sum_sim[sim_index] = cdp.q_pcs_norm_squared_sum_sim[sim_index] + cdp.q_factors_pcs_norm_squared_sum_sim[id][sim_index]**2
749 cdp.q_pcs_norm_squared_sum_sim[sim_index] = cdp.q_pcs_norm_squared_sum_sim[sim_index] / len(cdp.q_factors_pcs_norm_squared_sum_sim)
750 cdp.q_pcs_norm_squared_sum_sim[sim_index] = sqrt(cdp.q_pcs_norm_squared_sum_sim[sim_index])
751 else:
752 cdp.q_pcs_norm_squared_sum = 0.0
753 for id in cdp.q_factors_pcs_norm_squared_sum:
754 cdp.q_pcs_norm_squared_sum = cdp.q_pcs_norm_squared_sum + cdp.q_factors_pcs_norm_squared_sum[id]**2
755 cdp.q_pcs_norm_squared_sum = cdp.q_pcs_norm_squared_sum / len(cdp.q_factors_pcs_norm_squared_sum)
756 cdp.q_pcs_norm_squared_sum = sqrt(cdp.q_pcs_norm_squared_sum)
757
758
759 -def read(align_id=None, file=None, dir=None, file_data=None, spin_id_col=None, mol_name_col=None, res_num_col=None, res_name_col=None, spin_num_col=None, spin_name_col=None, data_col=None, error_col=None, sep=None, spin_id=None):
760 """Read the PCS data from file.
761
762 @param align_id: The alignment tensor ID string.
763 @type align_id: str
764 @param file: The name of the file to open.
765 @type file: str
766 @param dir: The directory containing the file (defaults to the current directory if None).
767 @type dir: str or None
768 @param file_data: An alternative opening a file, if the data already exists in the correct format. The format is a list of lists where the first index corresponds to the row and the second the column.
769 @type file_data: list of lists
770 @keyword spin_id_col: The column containing the spin ID strings. If supplied, the mol_name_col, res_name_col, res_num_col, spin_name_col, and spin_num_col arguments must be none.
771 @type spin_id_col: int or None
772 @keyword mol_name_col: The column containing the molecule name information. If supplied, spin_id_col must be None.
773 @type mol_name_col: int or None
774 @keyword res_name_col: The column containing the residue name information. If supplied, spin_id_col must be None.
775 @type res_name_col: int or None
776 @keyword res_num_col: The column containing the residue number information. If supplied, spin_id_col must be None.
777 @type res_num_col: int or None
778 @keyword spin_name_col: The column containing the spin name information. If supplied, spin_id_col must be None.
779 @type spin_name_col: int or None
780 @keyword spin_num_col: The column containing the spin number information. If supplied, spin_id_col must be None.
781 @type spin_num_col: int or None
782 @keyword data_col: The column containing the PCS data in Hz.
783 @type data_col: int or None
784 @keyword error_col: The column containing the PCS errors.
785 @type error_col: int or None
786 @keyword sep: The column separator which, if None, defaults to whitespace.
787 @type sep: str or None
788 @keyword spin_id: The spin ID string used to restrict data loading to a subset of all spins.
789 @type spin_id: None or str
790 """
791
792
793 check_pipe_setup(sequence=True)
794
795
796 if not exists_mol_res_spin_data():
797 raise RelaxNoSequenceError
798
799
800 if data_col == None and error_col == None:
801 raise RelaxError("One of either the data or error column must be supplied.")
802
803
804
805
806
807
808 mol_names = []
809 res_nums = []
810 res_names = []
811 spin_nums = []
812 spin_names = []
813 values = []
814 errors = []
815 for data in read_spin_data(file=file, dir=dir, file_data=file_data, spin_id_col=spin_id_col, mol_name_col=mol_name_col, res_num_col=res_num_col, res_name_col=res_name_col, spin_num_col=spin_num_col, spin_name_col=spin_name_col, data_col=data_col, error_col=error_col, sep=sep, spin_id=spin_id):
816
817 if data_col and error_col:
818 mol_name, res_num, res_name, spin_num, spin_name, value, error = data
819 elif data_col:
820 mol_name, res_num, res_name, spin_num, spin_name, value = data
821 error = None
822 else:
823 mol_name, res_num, res_name, spin_num, spin_name, error = data
824 value = None
825
826
827 if error == 0.0:
828 raise RelaxError("An invalid error value of zero has been encountered.")
829
830
831 id = generate_spin_id_unique(mol_name=mol_name, res_num=res_num, res_name=res_name, spin_num=spin_num, spin_name=spin_name)
832 spin = return_spin(spin_id=id)
833 if spin == None and spin_id and spin_id[0] == '@':
834 spin = return_spin(spin_id=id+spin_id)
835 if spin == None:
836 warn(RelaxNoSpinWarning(id))
837 continue
838
839
840 if data_col:
841
842 if not hasattr(spin, 'pcs'):
843 spin.pcs = {}
844
845
846 spin.pcs[align_id] = value
847
848
849 if error_col:
850
851 if not hasattr(spin, 'pcs_err'):
852 spin.pcs_err = {}
853
854
855 spin.pcs_err[align_id] = error
856
857
858 mol_names.append(mol_name)
859 res_nums.append(res_num)
860 res_names.append(res_name)
861 spin_nums.append(spin_num)
862 spin_names.append(spin_name)
863 values.append(value)
864 errors.append(error)
865
866
867 write_spin_data(file=sys.stdout, mol_names=mol_names, res_nums=res_nums, res_names=res_names, spin_nums=spin_nums, spin_names=spin_names, data=values, data_name='PCSs', error=errors, error_name='PCS_error')
868
869
870
871
872
873
874 if not len(values):
875 return
876
877
878 if not hasattr(cdp, 'align_ids'):
879 cdp.align_ids = []
880 if not hasattr(cdp, 'pcs_ids'):
881 cdp.pcs_ids = []
882
883
884 if align_id not in cdp.align_ids:
885 cdp.align_ids.append(align_id)
886 if align_id not in cdp.pcs_ids:
887 cdp.pcs_ids.append(align_id)
888
889
891 """Set up the data structures for optimisation using PCSs as base data sets.
892
893 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired.
894 @type sim_index: None or int
895 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity.
896 @type verbosity: int
897 @return: The assembled data structures for using PCSs as the base data for optimisation. These include:
898 - the PCS values.
899 - the unit vectors connecting the paramagnetic centre (the electron spin) to the spin.
900 - the PCS weight.
901 - the experimental temperatures.
902 - the spectrometer frequencies.
903 - pseudo_flags, the list of flags indicating if the interatomic data contains a pseudo-atom (as 1's and 0's).
904 @rtype: tuple of (numpy rank-2 float64 array, numpy rank-2 float64 array, numpy rank-2 float64 array, list of float, list of float, numpy rank-1 int32 array)
905 """
906
907
908 if verbosity:
909 print("\nPCS data counts:")
910
911
912 if not hasattr(cdp, 'paramagnetic_centre') and (hasattr(cdp, 'paramag_centre_fixed') and cdp.paramag_centre_fixed):
913 raise RelaxError("The paramagnetic centre has not yet been specified.")
914 if not hasattr(cdp, 'temperature'):
915 raise RelaxError("The experimental temperatures have not been set.")
916 if not hasattr(cdp, 'spectrometer_frq'):
917 raise RelaxError("The spectrometer frequencies of the experiments have not been set.")
918
919
920 setup_pseudoatom_pcs()
921
922
923 pcs = []
924 pcs_err = []
925 pcs_weight = []
926 temp = []
927 frq = []
928 pseudo_flags = []
929
930
931 for i in range(len(cdp.align_ids)):
932
933 align_id = cdp.align_ids[i]
934
935
936 if not opt_uses_align_data(align_id):
937 continue
938
939
940 pcs.append([])
941 pcs_err.append([])
942 pcs_weight.append([])
943
944
945 if align_id in cdp.temperature:
946 temp.append(cdp.temperature[align_id])
947
948
949 else:
950 raise RelaxError("The experimental temperature for the alignment ID '%s' has not been set." % align_id)
951
952
953 if align_id in cdp.spectrometer_frq:
954 frq.append(cdp.spectrometer_frq[align_id] * 2.0 * pi / periodic_table.gyromagnetic_ratio('1H'))
955
956
957 else:
958 raise RelaxError("The spectrometer frequency for the alignment ID '%s' has not been set." % align_id)
959
960
961 j = 0
962 for spin in spin_loop():
963
964 if not spin.select:
965 continue
966
967
968 if not hasattr(spin, 'pcs'):
969 continue
970
971
972 if align_id in spin.pcs:
973 if sim_index != None:
974 pcs[-1].append(spin.pcs_sim[align_id][sim_index])
975 else:
976 pcs[-1].append(spin.pcs[align_id])
977
978
979 if pcs[-1][-1] != None:
980 j += 1
981
982
983 else:
984 pcs[-1].append(None)
985
986
987 if hasattr(spin, 'pcs_err') and align_id in spin.pcs_err:
988 pcs_err[-1].append(spin.pcs_err[align_id])
989 else:
990 pcs_err[-1].append(None)
991
992
993 if hasattr(spin, 'pcs_weight') and align_id in spin.pcs_weight:
994 pcs_weight[-1].append(spin.pcs_weight[align_id])
995 else:
996 pcs_weight[-1].append(1.0)
997
998
999 if verbosity:
1000 print(" Alignment ID '%s': %i" % (align_id, j))
1001
1002
1003 for spin in spin_loop():
1004 if is_pseudoatom(spin):
1005 pseudo_flags.append(1)
1006 else:
1007 pseudo_flags.append(0)
1008
1009
1010 pcs = array(pcs, float64)
1011 pcs_err = array(pcs_err, float64)
1012 pcs_weight = array(pcs_weight, float64)
1013 pseudo_flags = array(pseudo_flags, int32)
1014
1015
1016 pcs = pcs * 1e-6
1017 pcs_err = pcs_err * 1e-6
1018
1019
1020 return pcs, pcs_err, pcs_weight, temp, frq, pseudo_flags
1021
1022
1023 -def set_errors(align_id=None, spin_id=None, sd=None):
1024 """Set the PCS errors if not already present.
1025
1026 @keyword align_id: The optional alignment tensor ID string.
1027 @type align_id: str
1028 @keyword spin_id: The optional spin ID string.
1029 @type spin_id: None or str
1030 @keyword sd: The PCS standard deviation in ppm.
1031 @type sd: float or int.
1032 """
1033
1034
1035 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True)
1036
1037
1038 if align_id:
1039 align_ids = [align_id]
1040 else:
1041 align_ids = cdp.pcs_ids
1042
1043
1044 for spin in spin_loop(spin_id):
1045
1046 if not spin.select:
1047 continue
1048
1049
1050 if not hasattr(spin, 'pcs') or (align_id and not align_id in spin.pcs):
1051 continue
1052
1053
1054 if not hasattr(spin, 'pcs_err'):
1055 spin.pcs_err = {}
1056
1057
1058 for id in align_ids:
1059 spin.pcs_err[id] = sd
1060
1061
1063 """Make sure that the spin systems are properly set up for pseudo-atoms and PCSs.
1064
1065 All spin data containers which are a member of a pseudo-atom will be deselected.
1066 """
1067
1068
1069 for pseudospin, pseudospin_id in spin_loop(return_id=True):
1070
1071 if not is_pseudoatom(pseudospin):
1072 return
1073
1074
1075 for spin, spin_id in pseudoatom_loop(pseudospin, return_id=True):
1076
1077 if spin.select:
1078 warn(RelaxWarning("Deselecting the '%s' spin as it is a member of the '%s' pseudo-atom system." % (spin_id, pseudospin_id)))
1079 spin.select = False
1080
1081
1082 -def structural_noise(align_id=None, rmsd=0.2, sim_num=1000, file=None, dir=None, force=False):
1083 """Determine the PCS error due to structural noise via simulation.
1084
1085 For the simulation the following must already be set up in the current data pipe:
1086
1087 - The position of the paramagnetic centre.
1088 - The alignment and magnetic susceptibility tensor.
1089
1090 The protocol for the simulation is as follows:
1091
1092 - The lanthanide or paramagnetic centre position will be fixed. Its motion is assumed to be on the femto- to pico- and nanosecond timescales. Hence the motion is averaged over the evolution of the PCS and can be ignored.
1093 - The positions of the nuclear spins will be randomised N times using a multivariate normal distribution.
1094 - The PCS for the randomised position will be back calculated.
1095 - The PCS standard deviation will be calculated from the N randomised PCS values.
1096
1097 The standard deviation will both be stored in the spin container data structure in the relax data store as well as being added to the already present PCS error (using variance addition). This will then be used in any optimisations involving the PCS.
1098
1099 If the alignment ID string is not supplied, the procedure will be applied to the PCS data from all alignments.
1100
1101
1102 @keyword align_id: The alignment tensor ID string.
1103 @type align_id: str
1104 @keyword rmsd: The atomic position RMSD, in Angstrom, to randomise the spin positions with for the simulations.
1105 @type rmsd: float
1106 @keyword sim_num: The number of simulations, N, to perform to determine the structural noise component of the PCS errors.
1107 @type sim_num: int
1108 @keyword file: The optional name of the Grace file to plot the structural errors verses the paramagnetic centre to spin distances.
1109 @type file: None or str
1110 @keyword dir: The directory name to place the Grace file into.
1111 @type dir: None or str
1112 @keyword force: A flag which if True will cause any pre-existing file to be overwritten.
1113 @type force: bool
1114 """
1115
1116
1117 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True, paramag_centre=True)
1118
1119
1120 if align_id:
1121 align_ids = [align_id]
1122 else:
1123 align_ids = cdp.align_ids
1124
1125
1126 grace_data = []
1127 unit_vect = zeros(3, float64)
1128 pcs = {}
1129 for id in align_ids:
1130 grace_data.append([])
1131 pcs[id] = zeros(sim_num, float64)
1132
1133
1134 cov = rmsd**2 * eye(3)
1135
1136
1137 print("Executing %i simulations for each spin system." % sim_num)
1138
1139
1140 for spin, spin_id in spin_loop(return_id=True):
1141
1142 if not spin.select:
1143 continue
1144
1145
1146 if not hasattr(spin, 'pcs'):
1147 continue
1148 if not hasattr(spin, 'pos'):
1149 continue
1150
1151
1152 print(spin_id)
1153
1154
1155 if type(spin.pos[0]) in [float, float64]:
1156 pos = spin.pos
1157 else:
1158 pos = zeros(3, float64)
1159 for i in range(len(spin.pos)):
1160 pos += spin.pos[i]
1161 pos = pos / len(spin.pos)
1162
1163
1164 orig_vect = pos - cdp.paramagnetic_centre
1165 orig_r = norm(orig_vect)
1166
1167
1168 new_pos = multivariate_normal(pos, cov, sim_num)
1169
1170
1171 for i in range(sim_num):
1172
1173 vect = new_pos[i] - cdp.paramagnetic_centre
1174 r = norm(vect)
1175 vect = vect / r
1176
1177
1178 for id in align_ids:
1179
1180 if id not in spin.pcs:
1181 continue
1182
1183
1184 dj = pcs_constant(cdp.temperature[id], cdp.spectrometer_frq[id] * 2.0 * pi / periodic_table.gyromagnetic_ratio('1H'), r/1e10)
1185
1186
1187 pcs[id][i] = pcs_tensor(dj, vect, cdp.align_tensors[get_tensor_index(id)].A) * 1e6
1188
1189
1190 if not hasattr(spin, 'pcs_struct_err'):
1191 spin.pcs_struct_err = {}
1192
1193
1194 align_index = 0
1195 for id in align_ids:
1196
1197 if id not in spin.pcs or spin.pcs[id] == None:
1198 align_index += 1
1199 continue
1200
1201
1202 sd = std(pcs[id])
1203
1204
1205 if id in spin.pcs_struct_err:
1206 warn(RelaxWarning("Removing the previous structural error value from the PCS error of the spin '%s' for the alignment ID '%s'." % (spin_id, id)))
1207 spin.pcs_err[id] = sqrt(spin.pcs_err[id]**2 - spin.pcs_struct_err[id]**2)
1208
1209
1210 spin.pcs_struct_err[id] = sd
1211
1212
1213 spin.pcs_err[id] = sqrt(spin.pcs_err[id]**2 + sd**2)
1214
1215
1216 grace_data[align_index].append([orig_r, sd, spin_id])
1217
1218
1219 align_index += 1
1220
1221
1222 if file:
1223
1224 file = open_write_file(file, dir, force)
1225
1226
1227 write_xy_header(format='grace', file=file, title="PCS structural noise", subtitle="%s Angstrom structural noise"%rmsd, data_type=['pcs_bc', 'pcs'], sets=[len(align_ids)], set_names=[align_ids], symbol_sizes=[[0.5]*len(align_ids)], linetype=[[0]*len(align_ids)], axis_labels=[["Ln\\S3+\\N to spin distance (Angstrom)", "PCS standard deviation (ppm)"]])
1228
1229
1230 write_xy_data(format='grace', data=[grace_data], file=file, graph_type='xy')
1231
1232
1233 -def weight(align_id=None, spin_id=None, weight=1.0):
1234 """Set optimisation weights on the PCS data.
1235
1236 @keyword align_id: The alignment tensor ID string.
1237 @type align_id: str
1238 @keyword spin_id: The spin ID string.
1239 @type spin_id: None or str
1240 @keyword weight: The optimisation weight. The higher the value, the more importance the PCS will have.
1241 @type weight: float or int.
1242 """
1243
1244
1245 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True)
1246
1247
1248 for spin in spin_loop(spin_id):
1249
1250 if not hasattr(spin, 'pcs_weight'):
1251 spin.pcs_weight = {}
1252
1253
1254 spin.pcs_weight[align_id] = weight
1255
1256
1257 -def write(align_id=None, file=None, dir=None, bc=False, force=False):
1258 """Display the PCS data corresponding to the alignment ID.
1259
1260 @keyword align_id: The alignment tensor ID string.
1261 @type align_id: str
1262 @keyword file: The file name or object to write to.
1263 @type file: str or file object
1264 @keyword dir: The name of the directory to place the file into (defaults to the current directory).
1265 @type dir: str
1266 @keyword bc: The back-calculation flag which if True will cause the back-calculated rather than measured data to be written.
1267 @type bc: bool
1268 @keyword force: A flag which if True will cause any pre-existing file to be overwritten.
1269 @type force: bool
1270 """
1271
1272
1273 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True)
1274
1275
1276 file = open_write_file(file, dir, force)
1277
1278
1279 mol_names = []
1280 res_nums = []
1281 res_names = []
1282 spin_nums = []
1283 spin_names = []
1284 values = []
1285 errors = []
1286 for spin, mol_name, res_num, res_name in spin_loop(full_info=True):
1287
1288 if not spin.select:
1289 continue
1290
1291
1292 if not bc and (not hasattr(spin, 'pcs') or not align_id in spin.pcs):
1293 continue
1294 elif bc and (not hasattr(spin, 'pcs_bc') or align_id not in spin.pcs_bc):
1295 continue
1296
1297
1298 mol_names.append(mol_name)
1299 res_nums.append(res_num)
1300 res_names.append(res_name)
1301 spin_nums.append(spin.num)
1302 spin_names.append(spin.name)
1303
1304
1305 if bc:
1306 values.append(spin.pcs_bc[align_id])
1307 else:
1308 values.append(spin.pcs[align_id])
1309
1310
1311 if hasattr(spin, 'pcs_err') and align_id in spin.pcs_err:
1312 errors.append(spin.pcs_err[align_id])
1313 else:
1314 errors.append(None)
1315
1316
1317 write_spin_data(file=file, mol_names=mol_names, res_nums=res_nums, res_names=res_names, spin_nums=spin_nums, spin_names=spin_names, data=values, data_name='PCSs', error=errors, error_name='PCS_error')
1318