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 functions for the handling of alignment tensors."""
24
25
26 from copy import deepcopy
27 from math import pi, sqrt
28 from numpy import arccos, complex128, dot, float64, inner, linalg, zeros
29 from numpy.linalg import norm
30 import sys
31 from warnings import warn
32
33
34 from data_store.align_tensor import AlignTensorList
35 from lib.alignment.alignment_tensor import calc_chi_tensor, kappa
36 from lib.errors import RelaxError, RelaxNoTensorError, RelaxTensorError, RelaxUnknownParamCombError, RelaxUnknownParamError
37 from lib.geometry.angles import wrap_angles
38 from lib.geometry.vectors import vector_angle_complex_conjugate
39 from lib.io import write_data
40 from lib.text.sectioning import section, subsection
41 from lib.warnings import RelaxWarning
42 import pipe_control
43 from pipe_control import pipes
44 from pipe_control.pipes import check_pipe
45
46
48 """Function for determining if alignment data exists in the current data pipe.
49
50 @keyword tensor: The alignment tensor ID string. If not supplied, then any alignment data will result in the function returning True.
51 @type tensor: str or None
52 @keyword pipe: The data pipe to search for data in. This defaults to the current data pipe if not supplied.
53 @type pipe: str or None
54 @return: The answer to the question.
55 @rtype: bool
56 """
57
58
59 if pipe == None:
60 pipe = pipes.cdp_name()
61
62
63 pipe = pipes.get_pipe(pipe)
64
65
66 if hasattr(pipe, 'align_tensors'):
67 if tensor == None:
68 return True
69 for data in pipe.align_tensors:
70 if data.name == tensor:
71 return True
72 else:
73 return False
74
75
77 """Determine if all alignment tensors are fixed.
78
79 @return: True if all tensors are fixed, False otherwise.
80 @rtype: bool
81 """
82
83
84 for i in range(len(cdp.align_tensors)):
85
86 if not cdp.align_tensors[i].fixed:
87 return False
88
89
90 return True
91
92
93 -def copy(tensor_from=None, pipe_from=None, tensor_to=None, pipe_to=None):
94 """Function for copying alignment tensor data from one data pipe to another.
95
96 @param tensor_from: The identification string of the alignment tensor to copy the data from.
97 @type tensor_from: str
98 @param pipe_from: The data pipe to copy the alignment tensor data from. This defaults to the current data pipe.
99 @type pipe_from: str
100 @param tensor_to: The identification string of the alignment tensor to copy the data to. If set to None, then the ID string will be set to the value of tensor_from.
101 @type tensor_to: str or None
102 @param pipe_to: The data pipe to copy the alignment tensor data to. This defaults to the current data pipe.
103 @type pipe_to: str
104 """
105
106
107 if tensor_from == tensor_to and pipe_from == None and pipe_to == None:
108 raise RelaxError("The pipe_from and pipe_to arguments cannot both be set to None when the tensor names are the same.")
109 elif pipe_from == None:
110 pipe_from = pipes.cdp_name()
111 elif pipe_to == None:
112 pipe_to = pipes.cdp_name()
113
114
115 if tensor_to == None:
116 tensor_to = tensor_from
117
118
119 check_pipe(pipe_from)
120 check_pipe(pipe_to)
121
122
123 dp_from = pipes.get_pipe(pipe_from)
124 dp_to = pipes.get_pipe(pipe_to)
125
126
127 if not align_data_exists(tensor_from, pipe_from):
128 raise RelaxNoTensorError('alignment')
129
130
131 if align_data_exists(tensor_to, pipe_to):
132 raise RelaxTensorError('alignment')
133
134
135 if not hasattr(dp_to, 'align_tensors'):
136 dp_to.align_tensors = AlignTensorList()
137
138
139 if tensor_to != None and tensor_to not in dp_to.align_tensors.names():
140 dp_to.align_tensors.add_item(tensor_to)
141
142
143 if tensor_from == None:
144
145 if tensor_to != tensor_from:
146 raise RelaxError("The tensor_to argument '%s' should not be supplied when tensor_from is None." % tensor_to)
147
148
149 align_ids = []
150
151
152 for i in range(len(dp_from.align_tensors)):
153 dp_to.align_tensors.append(deepcopy(dp_from.align_tensors[i]))
154 align_ids.append(dp_from.align_tensors[i].align_id)
155
156
157 else:
158
159 index_from = get_tensor_index(tensor=tensor_from, pipe=pipe_from)
160 index_to = get_tensor_index(tensor=tensor_to, pipe=pipe_to)
161
162
163 tensor = deepcopy(dp_from.align_tensors[index_from])
164 if index_to == None:
165 dp_to.align_tensors.append(tensor)
166 index_to = len(dp_to.align_tensors) - 1
167 else:
168 dp_to.align_tensors[index_to] = tensor
169
170
171 dp_to.align_tensors[index_to].set('name', tensor_to)
172
173
174 align_ids = [dp_from.align_tensors[index_from].align_id]
175
176
177 if not hasattr(dp_to, 'align_ids'):
178 dp_to.align_ids = []
179 for align_id in align_ids:
180 if align_id not in dp_to.align_ids:
181 dp_to.align_ids.append(align_id)
182
183
185 """Function for deleting alignment tensor data.
186
187 @param tensor: The alignment tensor identification string.
188 @type tensor: str or None
189 """
190
191
192 check_pipe()
193
194
195 if tensor and not align_data_exists(tensor):
196 raise RelaxNoTensorError('alignment')
197
198
199 if tensor:
200 tensors = [tensor]
201 else:
202 tensors = []
203 for i in range(len(cdp.align_tensors)):
204 tensors.append(cdp.align_tensors[i].name)
205
206
207 for tensor in tensors:
208
209 print("Removing the '%s' tensor." % tensor)
210
211
212 index = get_tensor_index(tensor=tensor)
213
214
215 cdp.align_tensors.pop(index)
216
217
218 if not len(cdp.align_tensors):
219 del(cdp.align_tensors)
220
221
223 """Function for displaying the alignment tensor.
224
225 @keyword tensor: The alignment tensor identification string.
226 @type tensor: str or None
227 """
228
229
230 check_pipe()
231
232
233 tensor_list = []
234 if not tensor:
235 for tensor_cont in cdp.align_tensors:
236 tensor_list.append(tensor_cont.name)
237 else:
238 tensor_list.append(tensor)
239
240
241 for tensor in tensor_list:
242
243 if not align_data_exists(tensor):
244 raise RelaxNoTensorError('alignment')
245
246
247 data = get_tensor_object(tensor)
248
249
250 section(file=sys.stdout, text="Tensor '%s'" % tensor, prespace=3, postspace=1)
251
252
253
254
255
256 subsection(file=sys.stdout, text="Saupe order matrix", prespace=0)
257
258
259 print("# 5D, rank-1 notation {Sxx, Syy, Sxy, Sxz, Syz}:")
260 print("[%25.12e, %25.12e, %25.12e, %25.12e, %25.12e]\n" % (data.Sxx, data.Syy, data.Sxy, data.Sxz, data.Syz))
261
262
263 print("# 5D, rank-1 notation {Szz, Sxx-yy, Sxy, Sxz, Syz} (the Pales default format).")
264 print("[%25.12e, %25.12e, %25.12e, %25.12e, %25.12e]\n" % (data.Szz, data.Sxxyy, data.Sxy, data.Sxz, data.Syz))
265
266
267 print("# 3D, rank-2 notation.")
268 print("%s" % (data.S))
269
270
271
272
273
274 subsection(file=sys.stdout, text="Alignment tensor", prespace=2)
275
276
277 print("# 5D, rank-1 notation {Axx, Ayy, Axy, Axz, Ayz}:")
278 print("[%25.12e, %25.12e, %25.12e, %25.12e, %25.12e]\n" % (data.Axx, data.Ayy, data.Axy, data.Axz, data.Ayz))
279
280
281 print("# 5D, rank-1 notation {Azz, Axx-yy, Axy, Axz, Ayz} (the Pales default format).")
282 print("[%25.12e, %25.12e, %25.12e, %25.12e, %25.12e]\n" % (data.Azz, data.Axxyy, data.Axy, data.Axz, data.Ayz))
283
284
285 print("# 3D, rank-2 notation.")
286 print("%s" % data.A)
287
288
289
290
291
292 subsection(file=sys.stdout, text="Probability tensor", prespace=2)
293
294
295 print("# 5D, rank-1 notation {Pxx, Pyy, Pxy, Pxz, Pyz}:")
296 print("[%25.12e, %25.12e, %25.12e, %25.12e, %25.12e]\n" % (data.Pxx, data.Pyy, data.Pxy, data.Pxz, data.Pyz))
297
298
299 print("# 5D, rank-1 notation {Pzz, Pxx-yy, Pxy, Pxz, Pyz}.")
300 print("[%25.12e, %25.12e, %25.12e, %25.12e, %25.12e]\n" % (data.Pzz, data.Pxxyy, data.Pxy, data.Pxz, data.Pyz))
301
302
303 print("# 3D, rank-2 notation.")
304 print("%s" % data.P)
305
306
307
308
309
310 subsection(file=sys.stdout, text="Magnetic susceptibility tensor", prespace=2)
311 chi_tensor = True
312
313
314 print("# The magnetic field strength (MHz):")
315 if hasattr(cdp, 'spectrometer_frq') and tensor in cdp.spectrometer_frq:
316 print("%s\n" % (cdp.spectrometer_frq[tensor] / 1e6))
317 else:
318 print("Not set.\n")
319 chi_tensor = False
320
321
322 print("# The temperature (K):")
323 if hasattr(cdp, 'temperature') and tensor in cdp.temperature:
324 print("%s\n" % cdp.temperature[tensor])
325 else:
326 print("Not set.\n")
327 chi_tensor = False
328
329
330 if not chi_tensor:
331 print("# The chi tensor:\nN/A.")
332
333
334 else:
335
336 chi_xx = calc_chi_tensor(data.Axx, cdp.spectrometer_frq[tensor], cdp.temperature[tensor])
337 chi_xy = calc_chi_tensor(data.Axy, cdp.spectrometer_frq[tensor], cdp.temperature[tensor])
338 chi_xz = calc_chi_tensor(data.Axz, cdp.spectrometer_frq[tensor], cdp.temperature[tensor])
339 chi_yy = calc_chi_tensor(data.Ayy, cdp.spectrometer_frq[tensor], cdp.temperature[tensor])
340 chi_yz = calc_chi_tensor(data.Ayz, cdp.spectrometer_frq[tensor], cdp.temperature[tensor])
341 chi_zz = calc_chi_tensor(data.Azz, cdp.spectrometer_frq[tensor], cdp.temperature[tensor])
342 chi_xxyy = calc_chi_tensor(data.Axxyy, cdp.spectrometer_frq[tensor], cdp.temperature[tensor])
343 chi = calc_chi_tensor(data.A, cdp.spectrometer_frq[tensor], cdp.temperature[tensor])
344
345
346 print("# 5D, rank-1 notation {chi_xx, chi_yy, chi_xy, chi_xz, chi_yz}:")
347 print("[%25.12e, %25.12e, %25.12e, %25.12e, %25.12e]\n" % (chi_xx, chi_yy, chi_xy, chi_xz, chi_yz))
348
349
350 print("# 5D, rank-1 notation {chi_zz, chi_xx-yy, chi_xy, chi_xz, chi_yz}.")
351 print("[%25.12e, %25.12e, %25.12e, %25.12e, %25.12e]\n" % (chi_zz, chi_xxyy, chi_xy, chi_xz, chi_yz))
352
353
354 print("# 3D, rank-2 notation.")
355 print("%s" % chi)
356
357
358
359
360
361 subsection(file=sys.stdout, text="Irreducible spherical tensor coefficients", prespace=2)
362
363
364 print("# The spherical harmonic decomposition weights are:")
365 print("# A0 = (4pi/5)^(1/2) Szz,")
366 print("# A+/-1 = +/- (8pi/15)^(1/2)(Sxz +/- iSyz),")
367 print("# A+/-2 = (2pi/15)^(1/2)(Sxx - Syy +/- 2iSxy).")
368
369
370 print("A-2 = %25.12e %25.12ei" % (data.Am2.real, data.Am2.imag))
371 print("A-1 = %25.12e %25.12ei" % (data.Am1.real, data.Am1.imag))
372 print("A0 = %25.12e" % data.A0)
373 print("A1 = %25.12e %25.12ei" % (data.A1.real, data.A1.imag))
374 print("A2 = %25.12e %25.12ei" % (data.A2.real, data.A2.imag))
375
376
377
378
379
380 subsection(file=sys.stdout, text="Eigensystem", prespace=2)
381
382
383 print("# Saupe order matrix eigenvalues {Sxx, Syy, Szz}.")
384 print("[%25.12e, %25.12e, %25.12e]\n" % (data.S_diag[0, 0], data.S_diag[1, 1], data.S_diag[2, 2]))
385 print("# Alignment tensor eigenvalues {Axx, Ayy, Azz}.")
386 print("[%25.12e, %25.12e, %25.12e]\n" % (data.A_diag[0, 0], data.A_diag[1, 1], data.A_diag[2, 2]))
387 print("# Probability tensor eigenvalues {Pxx, Pyy, Pzz}.")
388 print("[%25.12e, %25.12e, %25.12e]\n" % (data.P_diag[0, 0], data.P_diag[1, 1], data.P_diag[2, 2]))
389 if chi_tensor:
390 chi_diag = calc_chi_tensor(data.A_diag, cdp.spectrometer_frq[tensor], cdp.temperature[tensor])
391 print("# Magnetic susceptibility eigenvalues {chi_xx, chi_yy, chi_zz}.")
392 print("[%25.12e, %25.12e, %25.12e]\n" % (chi_diag[0, 0], chi_diag[1, 1], chi_diag[2, 2]))
393
394
395 print("# Eigenvector x.")
396 print("[%25.12f, %25.12f, %25.12f]\n" % (data.unit_x[0], data.unit_x[1], data.unit_x[2]))
397 print("# Eigenvector y.")
398 print("[%25.12f, %25.12f, %25.12f]\n" % (data.unit_y[0], data.unit_y[1], data.unit_y[2]))
399 print("# Eigenvector z.")
400 print("[%25.12f, %25.12f, %25.12f]\n" % (data.unit_z[0], data.unit_z[1], data.unit_z[2]))
401
402
403 print("# Rotation matrix.")
404 print("%s\n" % data.rotation)
405
406
407 print("# Euler angles in zyz notation {alpha, beta, gamma}.")
408 print("[%25.12f, %25.12f, %25.12f]\n" % (data.euler[0], data.euler[1], data.euler[2]))
409
410
411
412
413
414 subsection(file=sys.stdout, text="Geometric description", prespace=2)
415
416
417 print("# Generalized degree of order (GDO).")
418 print("GDO = %-25.12e\n" % gdo(data.A))
419
420
421 print("# Alignment tensor axial component (Aa = 3/2 * Azz, where Aii are the eigenvalues).")
422 print("Aa = %-25.12e\n" % data.Aa)
423
424
425 print("# Rhombic component (Ar = Axx - Ayy, where Aii are the eigenvalues).")
426 print("Ar = %-25.12e\n" % data.Ar)
427 print("# Rhombicity (R = Ar / Aa).")
428 print("R = %-25.12f\n" % data.R)
429 print("# Asymmetry parameter (eta = (Axx - Ayy) / Azz, where Aii are the eigenvalues).")
430 print("eta = %-25.12f\n" % data.eta)
431
432
433 if chi_tensor:
434
435 print("# Magnetic susceptibility axial parameter (chi_ax = chi_zz - (chi_xx + chi_yy)/2, where chi_ii are the eigenvalues).")
436 print("chi_ax = %-25.12e\n" % (chi_diag[2, 2] - (chi_diag[0, 0] + chi_diag[1, 1])/2.0))
437
438
439 print("# Magnetic susceptibility rhombicity parameter (chi_rh = chi_xx - chi_yy, where chi_ii are the eigenvalues).")
440 print("chi_rh = %-25.12e\n" % (chi_diag[0, 0] - chi_diag[1, 1]))
441
442
443 -def fix(id=None, fixed=True):
444 """Fix the alignment tensor during optimisation.
445
446 @keyword id: The alignment tensor ID string. If set to None, then all alignment tensors will be fixed.
447 @type id: str or None
448 @keyword fixed: If True, the alignment tensor will be fixed during optimisation. If False, the alignment tensors will be optimised.
449 @type fixed: bool
450 """
451
452
453 check_pipe()
454
455
456 for i in range(len(cdp.align_tensors)):
457
458 if id and cdp.align_tensors[i].name == id:
459 cdp.align_tensors[i].set_fixed(fixed)
460
461
462 if id == None:
463 cdp.align_tensors[i].set_fixed(fixed)
464
465
467 """Wrap the Euler angles and remove the glide reflection and translational symmetries.
468
469 Wrap the angles such that::
470
471 0 <= alpha <= 2pi,
472 0 <= beta <= pi,
473 0 <= gamma <= 2pi.
474
475 For the simulated values, the angles are wrapped as::
476
477 alpha - pi <= alpha_sim <= alpha + pi
478 beta - pi/2 <= beta_sim <= beta + pi/2
479 gamma - pi <= gamma_sim <= gamma + pi
480
481
482 @param sim_index: The simulation index. If set to None then the actual values will be folded
483 rather than the simulation values.
484 @type sim_index: int or None
485 """
486
487
488
489
490
491
492 alpha = cdp.align_tensors.alpha
493 beta = cdp.align_tensors.beta
494 gamma = cdp.align_tensors.gamma
495
496
497 if sim_index != None:
498 alpha_sim = cdp.align_tensors.alpha_sim[sim_index]
499 beta_sim = cdp.align_tensors.beta_sim[sim_index]
500 gamma_sim = cdp.align_tensors.gamma_sim[sim_index]
501
502
503 if sim_index == None:
504 cdp.align_tensors.set(param='alpha', value=wrap_angles(alpha, 0.0, 2.0*pi))
505 cdp.align_tensors.set(param='beta', value= wrap_angles(beta, 0.0, 2.0*pi))
506 cdp.align_tensors.set(param='gamma', value=wrap_angles(gamma, 0.0, 2.0*pi))
507
508
509 else:
510 cdp.align_tensors.set(param='alpha', value=wrap_angles(alpha_sim, alpha - pi, alpha + pi), category='sim', sim_index=sim_index)
511 cdp.align_tensors.set(param='beta', value= wrap_angles(beta_sim, beta - pi, beta + pi), category='sim', sim_index=sim_index)
512 cdp.align_tensors.set(param='gamma', value=wrap_angles(gamma_sim, gamma - pi, gamma + pi), category='sim', sim_index=sim_index)
513
514
515
516
517
518
519 if sim_index == None:
520
521 if cdp.align_tensors.beta >= pi:
522 cdp.align_tensors.set(param='alpha', value=pi - cdp.align_tensors.alpha)
523 cdp.align_tensors.set(param='beta', value=cdp.align_tensors.beta - pi)
524
525
526 else:
527
528 if cdp.align_tensors.beta_sim[sim_index] >= cdp.align_tensors.beta + pi/2.0:
529 cdp.align_tensors.set(param='alpha', value=pi - cdp.align_tensors.alpha_sim[sim_index], category='sim', sim_index=sim_index)
530 cdp.align_tensors.set(param='beta', value=cdp.align_tensors.beta_sim[sim_index] - pi, category='sim', sim_index=sim_index)
531 elif cdp.align_tensors.beta_sim[sim_index] <= cdp.align_tensors.beta - pi/2.0:
532 cdp.align_tensors.set(param='alpha', value=pi - cdp.align_tensors.alpha_sim[sim_index], category='sim', sim_index=sim_index)
533 cdp.align_tensors.set(param='beta', value=cdp.align_tensors.beta_sim[sim_index] + pi, category='sim', sim_index=sim_index)
534
535
537 """Calculate the generalized degree of order (GDO) for the given alignment tensor.
538
539 @param A: The alignment tensor.
540 @type A: rank-2, 3D numpy array
541 @return: The GDO value.
542 @rtype: float
543 """
544
545
546 gdo = sqrt(3.0/2.0) * norm(A)
547
548
549 return gdo
550
551
553 """Return the list of all alignment IDs.
554
555 @return: The list of all alignment IDs.
556 @rtype: list of str
557 """
558
559
560 if cdp == None:
561 return []
562
563
564 if not hasattr(cdp, 'align_ids'):
565 return []
566
567
568 return cdp.align_ids
569
570
572 """Return the list of all tensor IDs.
573
574 @return: The list of all tensor IDs.
575 @rtype: list of str
576 """
577
578
579 ids = []
580
581
582 if cdp == None:
583 return ids
584
585
586 if not hasattr(cdp, 'align_tensors'):
587 return ids
588
589
590 for i in range(len(cdp.align_tensors)):
591 if cdp.align_tensors[i].name != None:
592 ids.append(cdp.align_tensors[i].name)
593
594
595 return ids
596
597
599 """Function for returning the index corresponding to the 'tensor' argument.
600
601 @keyword tensor: The alignment tensor identification string.
602 @type tensor: str or None
603 @keyword align_id: Alternative to the tensor argument, used to return the tensor index for the tensors corresponding to the alignment ID string. If more than one tensor exists, then this will fail.
604 @type align_id: str or None
605 @keyword pipe: The data pipe to search for data in.
606 @type pipe: str
607 @return: The index corresponding to the 'tensor' arg.
608 @rtype: int
609 """
610
611
612 if pipe == None:
613 pipe = pipes.cdp_name()
614
615
616 dp = pipes.get_pipe(pipe)
617
618
619 index = None
620 count = 0
621
622
623 for i in range(len(dp.align_tensors)):
624
625 if tensor and dp.align_tensors[i].name == tensor:
626 index = i
627 count += 1
628
629
630 if align_id and hasattr(dp.align_tensors[i], 'align_id') and dp.align_tensors[i].align_id == align_id:
631 index = i
632 count += 1
633
634
635 if count == 0:
636 warn(RelaxWarning("No alignment tensors match the tensor name '%s' or alignment ID '%s' in the data pipe '%s'." % (tensor, align_id, pipe)))
637 return None
638
639
640 if count > 1:
641 warn(RelaxWarning("More than one alignment tensors matches the tensor name '%s' or alignment ID '%s' in the data pipe '%s'." % (tensor, align_id, pipe)))
642 return None
643
644
645 return index
646
647
649 """Return the AlignTensorData instance corresponding to the tensor ID.
650
651 @param tensor: The alignment tensor identification string.
652 @type tensor: str
653 @param pipe: The data pipe to search for data in.
654 @type pipe: str
655 @return: The alignment tensor object corresponding to the 'tensor' arg.
656 @rtype: AlignTensorData instance
657 """
658
659
660 if pipe == None:
661 pipe = pipes.cdp_name()
662
663
664 data = None
665
666
667 for i in range(len(cdp.align_tensors)):
668 if cdp.align_tensors[i].name == tensor:
669 data = cdp.align_tensors[i]
670
671
672 return data
673
674
676 """Return the unique AlignTensorData instance corresponding to the alignment ID.
677
678 @param align_id: The alignment ID for the unique tensor.
679 @type align_id: str
680 @return: The alignment tensor object corresponding to the 'tensor' arg.
681 @rtype: AlignTensorData instance
682 """
683
684
685 if pipe == None:
686 pipe = pipes.cdp_name()
687
688
689 data = None
690
691
692 count = 0
693 for i in range(len(cdp.align_tensors)):
694 if hasattr(cdp.align_tensors[i], 'align_id') and cdp.align_tensors[i].align_id == align_id:
695 data = cdp.align_tensors[i]
696 count += 1
697
698
699 if count > 1:
700 raise RelaxError("Multiple alignment tensors match the alignment ID '%s'." % align_id)
701
702 return data
703
704
705 -def init(tensor=None, align_id=None, params=None, scale=1.0, angle_units='deg', param_types=0, domain=None, errors=False):
706 """Function for initialising the alignment tensor.
707
708 @keyword tensor: The alignment tensor identification string.
709 @type tensor: str
710 @keyword align_id: The alignment ID string that the tensor corresponds to.
711 @type align_id: str or None
712 @keyword params: The alignment tensor parameters.
713 @type params: list of float or None
714 @keyword scale: The alignment tensor eigenvalue scaling value.
715 @type scale: float
716 @keyword angle_units: The units for the angle parameters (either 'deg' or 'rad').
717 @type angle_units: str
718 @keyword param_types: The type of parameters supplied. The flag values correspond to, 0: {Axx, Ayy, Axy, Axz, Ayz}, and 1: {Azz, Axx-yy, Axy, Axz, Ayz}.
719 @type param_types: int
720 @keyword domain: The domain label.
721 @type domain: str or None
722 @keyword errors: A flag which determines if the alignment tensor data or its errors are being input.
723 @type errors: bool
724 """
725
726
727 check_pipe()
728
729
730 if align_id == None:
731 raise RelaxError("The alignment ID must be given.")
732
733
734 valid_types = ['deg', 'rad']
735 if not angle_units in valid_types:
736 raise RelaxError("The alignment tensor 'angle_units' argument " + repr(angle_units) + " should be either 'deg' or 'rad'.")
737
738
739 if errors and (not hasattr(cdp, 'align_ids') or not align_id in cdp.align_ids):
740 raise RelaxNoTensorError('alignment')
741
742
743 if domain and (not hasattr(cdp, 'domain') or domain not in cdp.domain):
744 raise RelaxError("The domain '%s' has not been defined. Please use the domain user function." % domain)
745
746
747 if not hasattr(cdp, 'align_ids'):
748 cdp.align_ids = []
749 if align_id not in cdp.align_ids:
750 cdp.align_ids.append(align_id)
751
752
753 tensor_obj = None
754 if not errors:
755
756 if not hasattr(cdp, 'align_tensors'):
757 cdp.align_tensors = AlignTensorList()
758
759
760 if tensor == None or tensor not in cdp.align_tensors.names():
761 tensor_obj = cdp.align_tensors.add_item(tensor)
762
763
764 if not tensor_obj:
765 if tensor:
766 tensor_obj = get_tensor_object(tensor)
767 else:
768 tensor_obj = get_tensor_object_from_align(align_id)
769
770
771 if params:
772
773 if param_types == 0:
774
775 Sxx, Syy, Sxy, Sxz, Syz = params
776
777
778 Sxx = Sxx * scale
779 Syy = Syy * scale
780 Sxy = Sxy * scale
781 Sxz = Sxz * scale
782 Syz = Syz * scale
783
784
785 set(tensor=tensor_obj, value=[Sxx, Syy, Sxy, Sxz, Syz], param=['Sxx', 'Syy', 'Sxy', 'Sxz', 'Syz'], errors=errors)
786
787
788 elif param_types == 1:
789
790 Szz, Sxxyy, Sxy, Sxz, Syz = params
791
792
793 Szz = Szz * scale
794 Sxxyy = Sxxyy * scale
795 Sxy = Sxy * scale
796 Sxz = Sxz * scale
797 Syz = Syz * scale
798
799
800 set(tensor=tensor_obj, value=[Szz, Sxxyy, Sxy, Sxz, Syz], param=['Szz', 'Sxxyy', 'Sxy', 'Sxz', 'Syz'], errors=errors)
801
802
803 elif param_types == 2:
804
805 Axx, Ayy, Axy, Axz, Ayz = params
806
807
808 Axx = Axx * scale
809 Ayy = Ayy * scale
810 Axy = Axy * scale
811 Axz = Axz * scale
812 Ayz = Ayz * scale
813
814
815 set(tensor=tensor_obj, value=[Axx, Ayy, Axy, Axz, Ayz], param=['Axx', 'Ayy', 'Axy', 'Axz', 'Ayz'], errors=errors)
816
817
818 elif param_types == 3:
819
820 Azz, Axxyy, Axy, Axz, Ayz = params
821
822
823 Azz = Azz * scale
824 Axxyy = Axxyy * scale
825 Axy = Axy * scale
826 Axz = Axz * scale
827 Ayz = Ayz * scale
828
829
830 set(tensor=tensor_obj, value=[Azz, Axxyy, Axy, Axz, Ayz], param=['Azz', 'Axxyy', 'Axy', 'Axz', 'Ayz'], errors=errors)
831
832
833 elif param_types == 4:
834
835 Axx, Ayy, Axy, Axz, Ayz = params
836
837
838 r = None
839 for spin in spin_loop():
840
841 if r == None:
842 r = spin.r
843
844
845 if r != spin.r:
846 raise RelaxError("Not all spins have the same bond length.")
847
848
849 scale = scale / kappa() * r**3
850 Axx = Axx * scale
851 Ayy = Ayy * scale
852 Axy = Axy * scale
853 Axz = Axz * scale
854 Ayz = Ayz * scale
855
856
857 set(tensor=tensor_obj, value=[Axx, Ayy, Axy, Axz, Ayz], param=['Axx', 'Ayy', 'Axy', 'Axz', 'Ayz'], errors=errors)
858
859
860 elif param_types == 5:
861
862 Azz, Axxyy, Axy, Axz, Ayz = params
863
864
865 r = None
866 for spin in spin_loop():
867
868 if r == None:
869 r = spin.r
870
871
872 if r != spin.r:
873 raise RelaxError("Not all spins have the same bond length.")
874
875
876 scale = scale / kappa() * r**3
877 Azz = Azz * scale
878 Axxyy = Axxyy * scale
879 Axy = Axy * scale
880 Axz = Axz * scale
881 Ayz = Ayz * scale
882
883
884 set(tensor=tensor_obj, value=[Azz, Axxyy, Axy, Axz, Ayz], param=['Azz', 'Axxyy', 'Axy', 'Axz', 'Ayz'], errors=errors)
885
886
887 elif param_types == 6:
888
889 Pxx, Pyy, Pxy, Pxz, Pyz = params
890
891
892 Pxx = Pxx * scale
893 Pyy = Pyy * scale
894 Pxy = Pxy * scale
895 Pxz = Pxz * scale
896 Pyz = Pyz * scale
897
898
899 set(tensor=tensor_obj, value=[Pxx, Pyy, Pxy, Pxz, Pyz], param=['Pxx', 'Pyy', 'Pxy', 'Pxz', 'Pyz'], errors=errors)
900
901
902 elif param_types == 7:
903
904 Pzz, Pxxyy, Pxy, Pxz, Pyz = params
905
906
907 Pzz = Pzz * scale
908 Pxxyy = Pxxyy * scale
909 Pxy = Pxy * scale
910 Pxz = Pxz * scale
911 Pyz = Pyz * scale
912
913
914 set(tensor=tensor_obj, value=[Pzz, Pxxyy, Pxy, Pxz, Pyz], param=['Pzz', 'Pxxyy', 'Pxy', 'Pxz', 'Pyz'], errors=errors)
915
916
917 else:
918 raise RelaxUnknownParamCombError('param_types', param_types)
919
920
921 if domain:
922 set_domain(tensor=tensor, domain=domain)
923 if align_id:
924 tensor_obj.set(param='align_id', value=align_id)
925
926
927 -def matrix_angles(basis_set='matrix', tensors=None, angle_units='deg', precision=1):
928 """Function for calculating the inter-matrix angles between the alignment tensors.
929
930 The basis set defines how the angles are calculated:
931
932 - "matrix", the standard inter-matrix angle. The angle is calculated via the Euclidean inner product of the alignment matrices in rank-2, 3D form divided by the Frobenius norm ||A||_F of the matrices.
933 - "irreducible 5D", the irreducible 5D basis set {A-2, A-1, A0, A1, A2}.
934 - "unitary 5D", the unitary 5D basis set {Sxx, Syy, Sxy, Sxz, Syz}.
935 - "geometric 5D", the geometric 5D basis set {Szz, Sxxyy, Sxy, Sxz, Syz}. This is also the Pales standard notation.
936
937
938 @keyword basis_set: The basis set to use for calculating the inter-matrix angles. It can be one of "matrix", "irreducible 5D", "unitary 5D", or "geometric 5D".
939 @type basis_set: str
940 @keyword tensors: The list of alignment tensor IDs to calculate inter-matrix angles between. If None, all tensors will be used.
941 @type tensors: None or list of str
942 @keyword angle_units: The units for the angle parameters, either 'deg' or 'rad'.
943 @type angle_units: str
944 @keyword precision: The precision of the printed out angles. The number corresponds to the number of figures to print after the decimal point.
945 @type precision: int
946 """
947
948
949 allowed = ['matrix', 'unitary 9D', 'irreducible 5D', 'unitary 5D', 'geometric 5D']
950 if basis_set not in allowed:
951 raise RelaxError("The basis set of '%s' is not one of %s." % (basis_set, allowed))
952
953
954 if not hasattr(cdp, 'align_tensors') or len(cdp.align_tensors) == 0:
955 raise RelaxNoTensorError('alignment')
956
957
958 tensor_num = 0
959 for tensor in cdp.align_tensors:
960 if tensors and tensor.name not in tensors:
961 continue
962 tensor_num = tensor_num + 1
963
964
965 if basis_set == 'matrix':
966 matrix = zeros((tensor_num, 3, 3), float64)
967 elif basis_set == 'unitary 9D':
968 matrix = zeros((tensor_num, 9), float64)
969 elif basis_set in ['unitary 5D', 'geometric 5D']:
970 matrix = zeros((tensor_num, 5), float64)
971 elif basis_set in ['irreducible 5D']:
972 matrix = zeros((tensor_num, 5), complex128)
973 matrix_conj = zeros((tensor_num, 5), complex128)
974
975
976 i = 0
977 for tensor in cdp.align_tensors:
978
979 if tensors and tensor.name not in tensors:
980 continue
981
982
983 if basis_set == 'matrix':
984 matrix[i] = tensor.A
985
986
987 elif basis_set == 'unitary 9D':
988 matrix[i, 0] = tensor.Sxx
989 matrix[i, 1] = tensor.Sxy
990 matrix[i, 2] = tensor.Sxz
991 matrix[i, 3] = tensor.Sxy
992 matrix[i, 4] = tensor.Syy
993 matrix[i, 5] = tensor.Syz
994 matrix[i, 6] = tensor.Sxz
995 matrix[i, 7] = tensor.Syz
996 matrix[i, 8] = tensor.Szz
997
998
999 if basis_set == 'unitary 5D':
1000 matrix[i, 0] = tensor.Sxx
1001 matrix[i, 1] = tensor.Syy
1002 matrix[i, 2] = tensor.Sxy
1003 matrix[i, 3] = tensor.Sxz
1004 matrix[i, 4] = tensor.Syz
1005
1006
1007 if basis_set == 'irreducible 5D':
1008 matrix[i, 0] = tensor.Am2
1009 matrix[i, 1] = tensor.Am1
1010 matrix[i, 2] = tensor.A0
1011 matrix[i, 3] = tensor.A1
1012 matrix[i, 4] = tensor.A2
1013
1014
1015 matrix_conj[i, 0] = tensor.A2
1016 matrix_conj[i, 1] = -tensor.A1
1017 matrix_conj[i, 2] = tensor.A0
1018 matrix_conj[i, 3] = -tensor.Am1
1019 matrix_conj[i, 4] = tensor.Am2
1020
1021
1022 elif basis_set == 'geometric 5D':
1023 matrix[i, 0] = tensor.Szz
1024 matrix[i, 1] = tensor.Sxxyy
1025 matrix[i, 2] = tensor.Sxy
1026 matrix[i, 3] = tensor.Sxz
1027 matrix[i, 4] = tensor.Syz
1028
1029
1030 if basis_set in ['unitary 9D', 'unitary 5D', 'geometric 5D']:
1031 matrix[i] = matrix[i] / norm(matrix[i])
1032
1033
1034 i = i + 1
1035
1036
1037 cdp.align_tensors.angles = zeros((tensor_num, tensor_num), float64)
1038
1039
1040 if basis_set == 'matrix':
1041 sys.stdout.write("Standard inter-tensor matrix angles in degrees using the Euclidean inner product divided by the Frobenius norms (theta = arccos(<A1,A2>/(||A1||.||A2||)))")
1042 elif basis_set == 'irreducible 5D':
1043 sys.stdout.write("Inter-tensor vector angles in degrees for the irreducible 5D vectors {A-2, A-1, A0, A1, A2}")
1044 elif basis_set == 'unitary 9D':
1045 sys.stdout.write("Inter-tensor vector angles in degrees for the unitary 9D vectors {Sxx, Sxy, Sxz, Syx, Syy, Syz, Szx, Szy, Szz}")
1046 elif basis_set == 'unitary 5D':
1047 sys.stdout.write("Inter-tensor vector angles in degrees for the unitary 5D vectors {Sxx, Syy, Sxy, Sxz, Syz}")
1048 elif basis_set == 'geometric 5D':
1049 sys.stdout.write("Inter-tensor vector angles in degrees for the geometric 5D vectors {Szz, Sxx-yy, Sxy, Sxz, Syz}")
1050 sys.stdout.write(":\n\n")
1051
1052
1053 table = []
1054
1055
1056 table.append([''])
1057 for i in range(tensor_num):
1058
1059 if not tensors:
1060 if cdp.align_tensors[i].name == None:
1061 table[0].append(repr(i))
1062 else:
1063 table[0].append(cdp.align_tensors[i].name)
1064
1065
1066 else:
1067 table[0].append(tensors[i])
1068
1069
1070 for i in range(tensor_num):
1071
1072 if not tensors:
1073 if cdp.align_tensors[i].name == None:
1074 table.append([repr(i)])
1075 else:
1076 table.append([cdp.align_tensors[i].name])
1077 else:
1078 table.append([tensors[i]])
1079
1080
1081 for j in range(tensor_num):
1082
1083 if basis_set in ['unitary 9D', 'unitary 5D', 'geometric 5D']:
1084
1085 delta = dot(matrix[i], matrix[j])
1086
1087
1088 if delta > 1:
1089 delta = 1
1090
1091
1092 theta = arccos(delta)
1093
1094
1095 if basis_set in ['irreducible 5D']:
1096 theta = vector_angle_complex_conjugate(v1=matrix[i], v2=matrix[j], v1_conj=matrix_conj[i], v2_conj=matrix_conj[j])
1097
1098
1099 elif basis_set in ['matrix']:
1100
1101 nom = inner(matrix[i].flatten(), matrix[j].flatten())
1102
1103
1104 denom = norm(matrix[i]) * norm(matrix[j])
1105
1106
1107 ratio = nom / denom
1108 if ratio <= 1.0:
1109 theta = arccos(nom / denom)
1110 else:
1111 theta = 0.0
1112
1113
1114 cdp.align_tensors.angles[i, j] = theta
1115
1116
1117 angle = cdp.align_tensors.angles[i, j]
1118 if angle_units == 'deg':
1119 angle = angle * 180.0 / pi
1120 format = "%" + repr(7+precision) + "." + repr(precision) + "f"
1121 table[i+1].append(format % angle)
1122
1123
1124 write_data(out=sys.stdout, data=table)
1125
1126
1128 """Count the number of tensors.
1129
1130 @keyword skip_fixed: If set to True, then only the tensors without the fixed flag will be counted. If set to False, then all tensors will be counted.
1131 @type skip_fixed: bool
1132 @return: The number of tensors (excluding fixed tensors by default).
1133 @rtype: int
1134 """
1135
1136
1137 count = 0
1138
1139
1140 for tensor_cont in cdp.align_tensors:
1141
1142 if skip_fixed and tensor_cont.fixed:
1143 continue
1144
1145
1146 count += 1
1147
1148
1149 return count
1150
1151
1153 """Determine if the PCS or RDC data for the given alignment ID is needed for optimisation.
1154
1155 @keyword align_id: The optional alignment ID string.
1156 @type align_id: str
1157 @return: True if alignment data is to be used for optimisation, False otherwise.
1158 @rtype: bool
1159 """
1160
1161
1162 if not hasattr(cdp, 'align_ids'):
1163 return False
1164
1165
1166 if align_id:
1167 align_ids = [align_id]
1168 else:
1169 align_ids = cdp.align_ids
1170
1171
1172 for align_id in align_ids:
1173 if pipe_control.pcs.opt_uses_pcs(align_id) or pipe_control.rdc.opt_uses_rdc(align_id):
1174 return True
1175
1176
1177 return False
1178
1179
1181 """Determine if the given tensor is to be optimised.
1182
1183 @param tensor: The alignment tensor to check.
1184 @type tensor: AlignmentTensor object.
1185 @return: True if the tensor is to be optimised, False otherwise.
1186 @rtype: bool
1187 """
1188
1189
1190 ids = []
1191 if hasattr(cdp, 'rdc_ids'):
1192 ids += cdp.rdc_ids
1193 if hasattr(cdp, 'pcs_ids'):
1194 ids += cdp.pcs_ids
1195
1196
1197 if tensor.align_id not in ids:
1198 return False
1199
1200
1201 if tensor.fixed:
1202 return False
1203
1204
1205 return True
1206
1207
1208 -def reduction(full_tensor=None, red_tensor=None):
1209 """Specify which tensor is a reduction of which other tensor.
1210
1211 @param full_tensor: The full alignment tensor.
1212 @type full_tensor: str
1213 @param red_tensor: The reduced alignment tensor.
1214 @type red_tensor: str
1215 """
1216
1217
1218 match_full = False
1219 match_red = False
1220 i = 0
1221 for tensor_cont in cdp.align_tensors:
1222
1223 if tensor_cont.name == full_tensor:
1224 match_full = True
1225 index_full = i
1226 if tensor_cont.name == red_tensor:
1227 match_red = True
1228 index_red = i
1229
1230
1231 i = i + 1
1232
1233
1234 if not match_full:
1235 raise RelaxNoTensorError('alignment', full_tensor)
1236 if not match_red:
1237 raise RelaxNoTensorError('alignment', red_tensor)
1238
1239
1240 if not hasattr(cdp.align_tensors, 'reduction'):
1241 cdp.align_tensors.reduction = []
1242 cdp.align_tensors.reduction.append([index_full, index_red])
1243
1244
1246 """Return the tensor container for the given index, skipping fixed tensors if required.
1247
1248 @param index: The index of the tensor (if skip_fixed is True, then fixed tensors are not included in the index count).
1249 @type index: int
1250 @keyword skip_fixed: A flag which if True will exclude fixed tensors from the indexation.
1251 @type skip_fixed: bool
1252 @return: The tensor corresponding to the index.
1253 @rtype: data.align_tensor.AlignTensorData instance
1254 """
1255
1256
1257 count = 0
1258
1259
1260 for tensor_cont in cdp.align_tensors:
1261
1262 if skip_fixed and tensor_cont.fixed:
1263 continue
1264
1265
1266 if index == count:
1267 return tensor_cont
1268
1269
1270 count += 1
1271
1272
1273 return False
1274
1275
1276 -def set(tensor=None, value=None, param=None, errors=False):
1277 """Set the tensor.
1278
1279 @keyword tensor: The alignment tensor object.
1280 @type tensor: AlignTensorData instance
1281 @keyword value: The list of values to set the parameters to.
1282 @type value: list of float
1283 @keyword param: The list of parameter names.
1284 @type param: list of str
1285 @keyword errors: A flag which determines if the alignment tensor data or its errors are being
1286 input.
1287 @type errors: bool
1288 """
1289
1290
1291 geo_params = []
1292 geo_values = []
1293 orient_params = []
1294 orient_values = []
1295
1296
1297 for i in range(len(param)):
1298
1299 if param[i] == None:
1300 raise RelaxUnknownParamError("alignment tensor", 'None')
1301
1302
1303 if value[i] == None:
1304 value[i] = 0.0
1305
1306
1307 if param[i] in ['Sxx', 'Syy', 'Szz', 'Sxxyy', 'Sxy', 'Sxz', 'Syz', 'Axx', 'Ayy', 'Azz', 'Axxyy', 'Axy', 'Axz', 'Ayz', 'Pxx', 'Pyy', 'Pzz', 'Pxxyy', 'Pxy', 'Pxz', 'Pyz']:
1308 geo_params.append(param[i])
1309 geo_values.append(value[i])
1310
1311
1312 if param[i] in ['alpha', 'beta', 'gamma']:
1313 orient_params.append(param[i])
1314 orient_values.append(value[i])
1315
1316
1317
1318
1319
1320
1321 if len(geo_params) == 1:
1322
1323
1324
1325
1326 if geo_params[0] == 'Sxx':
1327 if errors:
1328 tensor.set(param='Sxx', value=geo_values[0], category='err')
1329 else:
1330 tensor.set(param='Sxx', value=geo_values[0])
1331
1332
1333 elif geo_params[0] == 'Syy':
1334 if errors:
1335 tensor.set(param='Syy', value=geo_values[0], category='err')
1336 else:
1337 tensor.set(param='Syy', value=geo_values[0])
1338
1339
1340 elif geo_params[0] == 'Sxy':
1341 if errors:
1342 tensor.set(param='Sxy', value=geo_values[0], category='err')
1343 else:
1344 tensor.set(param='Sxy', value=geo_values[0])
1345
1346
1347 elif geo_params[0] == 'Sxz':
1348 if errors:
1349 tensor.set(param='Sxz', value=geo_values[0], category='err')
1350 else:
1351 tensor.set(param='Sxz', value=geo_values[0])
1352
1353
1354 elif geo_params[0] == 'Syz':
1355 if errors:
1356 tensor.set(param='Syz', value=geo_values[0], category='err')
1357 else:
1358 tensor.set(param='Syz', value=geo_values[0])
1359
1360
1361
1362
1363
1364
1365 elif geo_params[0] == 'Axx':
1366 if errors:
1367 tensor.set(param='Sxx', value=3.0/2.0 * geo_values[0], category='err')
1368 else:
1369 tensor.set(param='Sxx', value=3.0/2.0 * geo_values[0])
1370
1371
1372 elif geo_params[0] == 'Ayy':
1373 if errors:
1374 tensor.set(param='Syy', value=3.0/2.0 * geo_values[0], category='err')
1375 else:
1376 tensor.set(param='Syy', value=3.0/2.0 * geo_values[0])
1377
1378
1379 elif geo_params[0] == 'Axy':
1380 if errors:
1381 tensor.set(param='Sxy', value=3.0/2.0 * geo_values[0], category='err')
1382 else:
1383 tensor.set(param='Sxy', value=3.0/2.0 * geo_values[0])
1384
1385
1386 elif geo_params[0] == 'Axz':
1387 if errors:
1388 tensor.set(param='Sxz', value=3.0/2.0 * geo_values[0], category='err')
1389 else:
1390 tensor.set(param='Sxz', value=3.0/2.0 * geo_values[0])
1391
1392
1393 elif geo_params[0] == 'Ayz':
1394 if errors:
1395 tensor.set(param='Syz', value=3.0/2.0 * geo_values[0], category='err')
1396 else:
1397 tensor.set(param='Syz', value=3.0/2.0 * geo_values[0])
1398
1399
1400
1401
1402
1403
1404 elif geo_params[0] == 'Pxx':
1405 if errors:
1406 tensor.set(param='Sxx', value=3.0/2.0 * (geo_values[0] - 1.0/3.0), category='err')
1407 else:
1408 tensor.set(param='Sxx', value=3.0/2.0 * (geo_values[0] - 1.0/3.0))
1409
1410
1411 elif geo_params[0] == 'Pyy':
1412 if errors:
1413 tensor.set(param='Syy', value=3.0/2.0 * (geo_values[0] - 1.0/3.0), category='err')
1414 else:
1415 tensor.set(param='Syy', value=3.0/2.0 * (geo_values[0] - 1.0/3.0))
1416
1417
1418 elif geo_params[0] == 'Pxy':
1419 if errors:
1420 tensor.set(param='Sxy', value=3.0/2.0 * geo_values[0], category='err')
1421 else:
1422 tensor.set(param='Sxy', value=3.0/2.0 * geo_values[0])
1423
1424
1425 elif geo_params[0] == 'Pxz':
1426 if errors:
1427 tensor.set(param='Sxz', value=3.0/2.0 * geo_values[0], category='err')
1428 else:
1429 tensor.set(param='Sxz', value=3.0/2.0 * geo_values[0])
1430
1431
1432 elif geo_params[0] == 'Pyz':
1433 if errors:
1434 tensor.set(param='Syz', value=3.0/2.0 * geo_values[0], category='err')
1435 else:
1436 tensor.set(param='Syz', value=3.0/2.0 * geo_values[0])
1437
1438
1439 else:
1440 raise RelaxError("The geometric alignment parameter " + repr(geo_params[0]) + " cannot be set.")
1441
1442
1443 elif len(geo_params) == 5:
1444
1445 if geo_params.count('Sxx') == 1 and geo_params.count('Syy') == 1 and geo_params.count('Sxy') == 1 and geo_params.count('Sxz') == 1 and geo_params.count('Syz') == 1:
1446
1447 Sxx = geo_values[geo_params.index('Sxx')]
1448 Syy = geo_values[geo_params.index('Syy')]
1449 Sxy = geo_values[geo_params.index('Sxy')]
1450 Sxz = geo_values[geo_params.index('Sxz')]
1451 Syz = geo_values[geo_params.index('Syz')]
1452
1453
1454 if errors:
1455 tensor.set(param='Axx', value=2.0/3.0 * Sxx, category='err')
1456 tensor.set(param='Ayy', value=2.0/3.0 * Syy, category='err')
1457 tensor.set(param='Axy', value=2.0/3.0 * Sxy, category='err')
1458 tensor.set(param='Axz', value=2.0/3.0 * Sxz, category='err')
1459 tensor.set(param='Ayz', value=2.0/3.0 * Syz, category='err')
1460 else:
1461 tensor.set(param='Axx', value=2.0/3.0 * Sxx)
1462 tensor.set(param='Ayy', value=2.0/3.0 * Syy)
1463 tensor.set(param='Axy', value=2.0/3.0 * Sxy)
1464 tensor.set(param='Axz', value=2.0/3.0 * Sxz)
1465 tensor.set(param='Ayz', value=2.0/3.0 * Syz)
1466
1467
1468 elif geo_params.count('Szz') == 1 and geo_params.count('Sxxyy') == 1 and geo_params.count('Sxy') == 1 and geo_params.count('Sxz') == 1 and geo_params.count('Syz') == 1:
1469
1470 Szz = geo_values[geo_params.index('Szz')]
1471 Sxxyy = geo_values[geo_params.index('Sxxyy')]
1472 Sxy = geo_values[geo_params.index('Sxy')]
1473 Sxz = geo_values[geo_params.index('Sxz')]
1474 Syz = geo_values[geo_params.index('Syz')]
1475
1476
1477 if errors:
1478 tensor.set(param='Axx', value=2.0/3.0 * -0.5*(Szz-Sxxyy), category='err')
1479 tensor.set(param='Ayy', value=2.0/3.0 * -0.5*(Szz+Sxxyy), category='err')
1480 tensor.set(param='Axy', value=2.0/3.0 * Sxy, category='err')
1481 tensor.set(param='Axz', value=2.0/3.0 * Sxz, category='err')
1482 tensor.set(param='Ayz', value=2.0/3.0 * Syz, category='err')
1483 else:
1484 tensor.set(param='Axx', value=2.0/3.0 * -0.5*(Szz-Sxxyy))
1485 tensor.set(param='Ayy', value=2.0/3.0 * -0.5*(Szz+Sxxyy))
1486 tensor.set(param='Axy', value=2.0/3.0 * Sxy)
1487 tensor.set(param='Axz', value=2.0/3.0 * Sxz)
1488 tensor.set(param='Ayz', value=2.0/3.0 * Syz)
1489
1490
1491 elif geo_params.count('Axx') == 1 and geo_params.count('Ayy') == 1 and geo_params.count('Axy') == 1 and geo_params.count('Axz') == 1 and geo_params.count('Ayz') == 1:
1492
1493 Axx = geo_values[geo_params.index('Axx')]
1494 Ayy = geo_values[geo_params.index('Ayy')]
1495 Axy = geo_values[geo_params.index('Axy')]
1496 Axz = geo_values[geo_params.index('Axz')]
1497 Ayz = geo_values[geo_params.index('Ayz')]
1498
1499
1500 if errors:
1501 tensor.set(param='Axx', value=Axx, category='err')
1502 tensor.set(param='Ayy', value=Ayy, category='err')
1503 tensor.set(param='Axy', value=Axy, category='err')
1504 tensor.set(param='Axz', value=Axz, category='err')
1505 tensor.set(param='Ayz', value=Ayz, category='err')
1506 else:
1507 tensor.set(param='Axx', value=Axx)
1508 tensor.set(param='Ayy', value=Ayy)
1509 tensor.set(param='Axy', value=Axy)
1510 tensor.set(param='Axz', value=Axz)
1511 tensor.set(param='Ayz', value=Ayz)
1512
1513
1514 elif geo_params.count('Azz') == 1 and geo_params.count('Axxyy') == 1 and geo_params.count('Axy') == 1 and geo_params.count('Axz') == 1 and geo_params.count('Ayz') == 1:
1515
1516 Azz = geo_values[geo_params.index('Azz')]
1517 Axxyy = geo_values[geo_params.index('Axxyy')]
1518 Axy = geo_values[geo_params.index('Axy')]
1519 Axz = geo_values[geo_params.index('Axz')]
1520 Ayz = geo_values[geo_params.index('Ayz')]
1521
1522
1523 if errors:
1524 tensor.set(param='Axx', value=-0.5*(Azz-Axxyy), category='err')
1525 tensor.set(param='Ayy', value=-0.5*(Azz+Axxyy), category='err')
1526 tensor.set(param='Axy', value=Axy, category='err')
1527 tensor.set(param='Axz', value=Axz, category='err')
1528 tensor.set(param='Ayz', value=Ayz, category='err')
1529 else:
1530 tensor.set(param='Axx', value=-0.5*(Azz-Axxyy))
1531 tensor.set(param='Ayy', value=-0.5*(Azz+Axxyy))
1532 tensor.set(param='Axy', value=Axy)
1533 tensor.set(param='Axz', value=Axz)
1534 tensor.set(param='Ayz', value=Ayz)
1535
1536
1537 elif geo_params.count('Pxx') == 1 and geo_params.count('Pyy') == 1 and geo_params.count('Pxy') == 1 and geo_params.count('Pxz') == 1 and geo_params.count('Pyz') == 1:
1538
1539 Pxx = geo_values[geo_params.index('Pxx')]
1540 Pyy = geo_values[geo_params.index('Pyy')]
1541 Pxy = geo_values[geo_params.index('Pxy')]
1542 Pxz = geo_values[geo_params.index('Pxz')]
1543 Pyz = geo_values[geo_params.index('Pyz')]
1544
1545
1546 if errors:
1547 tensor.set(param='Axx', value=Pxx - 1.0/3.0, category='err')
1548 tensor.set(param='Ayy', value=Pyy - 1.0/3.0, category='err')
1549 tensor.set(param='Axy', value=Pxy, category='err')
1550 tensor.set(param='Axz', value=Pxz, category='err')
1551 tensor.set(param='Ayz', value=Pyz, category='err')
1552 else:
1553 tensor.set(param='Axx', value=Pxx - 1.0/3.0)
1554 tensor.set(param='Ayy', value=Pyy - 1.0/3.0)
1555 tensor.set(param='Axy', value=Pxy)
1556 tensor.set(param='Axz', value=Pxz)
1557 tensor.set(param='Ayz', value=Pyz)
1558
1559
1560 elif geo_params.count('Pzz') == 1 and geo_params.count('Pxxyy') == 1 and geo_params.count('Pxy') == 1 and geo_params.count('Pxz') == 1 and geo_params.count('Pyz') == 1:
1561
1562 Pzz = geo_values[geo_params.index('Pzz')]
1563 Pxxyy = geo_values[geo_params.index('Pxxyy')]
1564 Pxy = geo_values[geo_params.index('Pxy')]
1565 Pxz = geo_values[geo_params.index('Pxz')]
1566 Pyz = geo_values[geo_params.index('Pyz')]
1567
1568
1569 if errors:
1570 tensor.set(param='Axx', value=-0.5*(Pzz-Pxxyy) - 1.0/3.0, category='err')
1571 tensor.set(param='Ayy', value=-0.5*(Pzz+Pxxyy) - 1.0/3.0, category='err')
1572 tensor.set(param='Axy', value=Pxy, category='err')
1573 tensor.set(param='Axz', value=Pxz, category='err')
1574 tensor.set(param='Ayz', value=Pyz, category='err')
1575 else:
1576 tensor.set(param='Axx', value=-0.5*(Pzz-Pxxyy) - 1.0/3.0)
1577 tensor.set(param='Ayy', value=-0.5*(Pzz+Pxxyy) - 1.0/3.0)
1578 tensor.set(param='Axy', value=Pxy)
1579 tensor.set(param='Axz', value=Pxz)
1580 tensor.set(param='Ayz', value=Pyz)
1581
1582
1583 else:
1584 raise RelaxUnknownParamCombError('geometric parameter set', geo_params)
1585
1586
1587
1588 else:
1589 raise RelaxUnknownParamCombError('geometric parameter set', geo_params)
1590
1591
1592
1593
1594
1595
1596 if len(orient_params) == 1:
1597
1598 if orient_params[0] == 'alpha':
1599 if errors:
1600 tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')], category='err')
1601 else:
1602 tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')])
1603
1604
1605 elif orient_params[0] == 'beta':
1606 if errors:
1607 tensor.set(param='beta', value=orient_values[orient_params.index('beta')], category='err')
1608 else:
1609 tensor.set(param='beta', value=orient_values[orient_params.index('beta')])
1610
1611
1612 elif orient_params[0] == 'gamma':
1613 if errors:
1614 tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')], category='err')
1615 else:
1616 tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')])
1617
1618
1619 elif len(orient_params) == 2:
1620
1621 if orient_params.count('alpha') == 1 and orient_params.count('beta') == 1:
1622 if errors:
1623 tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')], category='err')
1624 tensor.set(param='beta', value=orient_values[orient_params.index('beta')], category='err')
1625 else:
1626 tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')])
1627 tensor.set(param='beta', value=orient_values[orient_params.index('beta')])
1628
1629
1630 if orient_params.count('alpha') == 1 and orient_params.count('gamma') == 1:
1631 if errors:
1632 tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')], category='err')
1633 tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')], category='err')
1634 else:
1635 tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')])
1636 tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')])
1637
1638
1639 if orient_params.count('beta') == 1 and orient_params.count('gamma') == 1:
1640 if errors:
1641 tensor.set(param='beta', value=orient_values[orient_params.index('beta')], category='err')
1642 tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')], category='err')
1643 else:
1644 tensor.set(param='beta', value=orient_values[orient_params.index('beta')])
1645 tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')])
1646
1647
1648 else:
1649 raise RelaxUnknownParamCombError('orientational parameter set', orient_params)
1650
1651
1652 elif len(orient_params) == 3:
1653
1654 if orient_params.count('alpha') == 1 and orient_params.count('beta') == 1:
1655 if errors:
1656 tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')], category='err')
1657 tensor.set(param='beta', value=orient_values[orient_params.index('beta')], category='err')
1658 tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')], category='err')
1659 else:
1660 tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')])
1661 tensor.set(param='beta', value=orient_values[orient_params.index('beta')])
1662 tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')])
1663
1664
1665 else:
1666 raise RelaxUnknownParamCombError('orientational parameter set', orient_params)
1667
1668
1669 elif len(orient_params) > 3:
1670 raise RelaxUnknownParamCombError('orientational parameter set', orient_params)
1671
1672
1673
1674
1675
1676 if orient_params:
1677 fold_angles()
1678
1679
1681 """Set the align ID string for the given tensor.
1682
1683 @keyword tensor: The alignment tensor label.
1684 @type tensor: str
1685 @keyword align_id: The alignment ID string.
1686 @type align_id: str
1687 """
1688
1689
1690 match = False
1691 for tensor_cont in cdp.align_tensors:
1692
1693 if tensor_cont.name == tensor:
1694 tensor_cont.align_id = align_id
1695 match = True
1696
1697
1698 if not match:
1699 raise RelaxNoTensorError('alignment', tensor)
1700
1701
1702 -def set_domain(tensor=None, domain=None):
1703 """Set the domain label for the given tensor.
1704
1705 @param tensor: The alignment tensor label.
1706 @type tensor: str
1707 @param domain: The domain label.
1708 @type domain: str
1709 """
1710
1711
1712 if not hasattr(cdp, 'domain') or domain not in cdp.domain:
1713 raise RelaxError("The domain '%s' has not been defined. Please use the domain user function." % domain)
1714
1715
1716 match = False
1717 for tensor_cont in cdp.align_tensors:
1718
1719 if tensor_cont.name == tensor:
1720 tensor_cont.set(param='domain', value=domain)
1721 match = True
1722
1723
1724 if not match:
1725 raise RelaxNoTensorError('alignment', tensor)
1726
1727
1728 -def svd(basis_set='irreducible 5D', tensors=None, precision=1):
1729 r"""Calculate the singular values of all the loaded tensors.
1730
1731 The basis set can be set to one of:
1732
1733 - 'irreducible 5D', the irreducible 5D basis set {A-2, A-1, A0, A1, A2}. This is a linear map, hence angles are preserved.
1734 - 'unitary 9D', the unitary 9D basis set {Sxx, Sxy, Sxz, Syx, Syy, Syz, Szx, Szy, Szz}. This is a linear map, hence angles are preserved.
1735 - 'unitary 5D', the unitary 5D basis set {Sxx, Syy, Sxy, Sxz, Syz}. This is a non-linear map, hence angles are not preserved.
1736 - 'geometric 5D', the geometric 5D basis set {Szz, Sxxyy, Sxy, Sxz, Syz}. This is a non-linear map, hence angles are not preserved. This is also the Pales standard notation.
1737
1738 If the selected basis set is the default of 'irreducible 5D', the matrix on which SVD will be performed will be::
1739
1740 | A-2(1) A-1(1) A0(1) A1(1) A2(1) |
1741 | A-2(2) A-1(2) A0(2) A1(2) A2(2) |
1742 | A-2(3) A-1(3) A0(3) A1(3) A2(3) |
1743 | . . . . . |
1744 | . . . . . |
1745 | . . . . . |
1746 | A-2(N) A-1(N) A0(N) A1(N) A2(N) |
1747
1748 If the selected basis set is 'unitary 9D', the matrix on which SVD will be performed will be::
1749
1750 | Sxx1 Sxy1 Sxz1 Syx1 Syy1 Syz1 Szx1 Szy1 Szz1 |
1751 | Sxx2 Sxy2 Sxz2 Syx2 Syy2 Syz2 Szx2 Szy2 Szz2 |
1752 | Sxx3 Sxy3 Sxz3 Syx3 Syy3 Syz3 Szx3 Szy3 Szz3 |
1753 | . . . . . . . . . |
1754 | . . . . . . . . . |
1755 | . . . . . . . . . |
1756 | SxxN SxyN SxzN SyxN SyyN SyzN SzxN SzyN SzzN |
1757
1758 Otherwise if the selected basis set is 'unitary 5D', the matrix for SVD is::
1759
1760 | Sxx1 Syy1 Sxy1 Sxz1 Syz1 |
1761 | Sxx2 Syy2 Sxy2 Sxz2 Syz2 |
1762 | Sxx3 Syy3 Sxy3 Sxz3 Syz3 |
1763 | . . . . . |
1764 | . . . . . |
1765 | . . . . . |
1766 | SxxN SyyN SxyN SxzN SyzN |
1767
1768 Or if the selected basis set is 'geometric 5D', the stretching and skewing parameters Szz and Sxx-yy will be used instead and the matrix is::
1769
1770 | Szz1 Sxxyy1 Sxy1 Sxz1 Syz1 |
1771 | Szz2 Sxxyy2 Sxy2 Sxz2 Syz2 |
1772 | Szz3 Sxxyy3 Sxy3 Sxz3 Syz3 |
1773 | . . . . . |
1774 | . . . . . |
1775 | . . . . . |
1776 | SzzN SxxyyN SxyN SxzN SyzN |
1777
1778 For the irreducible basis set, the Am components are defined as::
1779
1780 / 4pi \ 1/2
1781 A0 = | --- | Szz ,
1782 \ 5 /
1783
1784 / 8pi \ 1/2
1785 A+/-1 = +/- | --- | (Sxz +/- iSyz) ,
1786 \ 15 /
1787
1788 / 2pi \ 1/2
1789 A+/-2 = | --- | (Sxx - Syy +/- 2iSxy) .
1790 \ 15 /
1791
1792 The relationships between the geometric and unitary basis sets are::
1793
1794 Szz = - Sxx - Syy,
1795 Sxxyy = Sxx - Syy,
1796
1797 The SVD values and condition number are dependant upon the basis set chosen.
1798
1799
1800 @keyword basis_set: The basis set to use for the SVD. This can be one of "irreducible 5D", "unitary 9D", "unitary 5D" or "geometric 5D".
1801 @type basis_set: str
1802 @keyword tensors: The list of alignment tensor IDs to calculate inter-matrix angles between. If None, all tensors will be used.
1803 @type tensors: None or list of str
1804 @keyword precision: The precision of the printed out angles. The number corresponds to the number of figures to print after the decimal point.
1805 @type precision: int
1806 """
1807
1808
1809 allowed = ['irreducible 5D', 'unitary 9D', 'unitary 5D', 'geometric 5D']
1810 if basis_set not in allowed:
1811 raise RelaxError("The basis set of '%s' is not one of %s." % (basis_set, allowed))
1812
1813
1814 if not hasattr(cdp, 'align_tensors') or len(cdp.align_tensors) == 0:
1815 raise RelaxNoTensorError('alignment')
1816
1817
1818 tensor_num = 0
1819 for tensor in cdp.align_tensors:
1820 if tensors and tensor.name not in tensors:
1821 continue
1822 tensor_num = tensor_num + 1
1823
1824
1825 if basis_set in ['unitary 9D']:
1826 matrix = zeros((tensor_num, 9), float64)
1827 elif basis_set in ['irreducible 5D']:
1828 matrix = zeros((tensor_num, 5), complex128)
1829 else:
1830 matrix = zeros((tensor_num, 5), float64)
1831
1832
1833 i = 0
1834 for tensor in cdp.align_tensors:
1835
1836 if tensors and tensor.name not in tensors:
1837 continue
1838
1839
1840 if basis_set == 'irreducible 5D':
1841 matrix[i, 0] = tensor.Am2
1842 matrix[i, 1] = tensor.Am1
1843 matrix[i, 2] = tensor.A0
1844 matrix[i, 3] = tensor.A1
1845 matrix[i, 4] = tensor.A2
1846
1847
1848 elif basis_set == 'unitary 9D':
1849 matrix[i, 0] = tensor.Sxx
1850 matrix[i, 1] = tensor.Sxy
1851 matrix[i, 2] = tensor.Sxz
1852 matrix[i, 3] = tensor.Sxy
1853 matrix[i, 4] = tensor.Syy
1854 matrix[i, 5] = tensor.Syz
1855 matrix[i, 6] = tensor.Sxz
1856 matrix[i, 7] = tensor.Syz
1857 matrix[i, 8] = tensor.Szz
1858
1859
1860 elif basis_set == 'unitary 5D':
1861 matrix[i, 0] = tensor.Sxx
1862 matrix[i, 1] = tensor.Syy
1863 matrix[i, 2] = tensor.Sxy
1864 matrix[i, 3] = tensor.Sxz
1865 matrix[i, 4] = tensor.Syz
1866
1867
1868 elif basis_set == 'geometric 5D':
1869 matrix[i, 0] = tensor.Szz
1870 matrix[i, 1] = tensor.Sxxyy
1871 matrix[i, 2] = tensor.Sxy
1872 matrix[i, 3] = tensor.Sxz
1873 matrix[i, 4] = tensor.Syz
1874
1875
1876 i = i + 1
1877
1878
1879 u, s, vh = linalg.svd(matrix)
1880
1881
1882 cdp.align_tensors.singular_vals = s
1883
1884
1885 cdp.align_tensors.cond_num = s[0] / s[-1]
1886
1887
1888 if basis_set == 'irreducible 5D':
1889 sys.stdout.write("SVD for the irreducible 5D vectors {A-2, A-1, A0, A1, A2}.\n")
1890 elif basis_set == 'unitary 9D':
1891 sys.stdout.write("SVD for the unitary 9D vectors {Sxx, Sxy, Sxz, Syx, Syy, Syz, Szx, Szy, Szz}.\n")
1892 elif basis_set == 'unitary 5D':
1893 sys.stdout.write("SVD for the unitary 5D vectors {Sxx, Syy, Sxy, Sxz, Syz}.\n")
1894 elif basis_set == 'geometric 5D':
1895 sys.stdout.write("SVD for the geometric 5D vectors {Szz, Sxx-yy, Sxy, Sxz, Syz}.\n")
1896 sys.stdout.write("\nSingular values:\n")
1897 for val in s:
1898 format = " %." + repr(precision) + "e\n"
1899 sys.stdout.write(format % val)
1900 format = "\nCondition number: %." + repr(precision) + "f\n"
1901 sys.stdout.write(format % cdp.align_tensors.cond_num)
1902