1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 """Module for the handling of Frame Order."""
24
25
26 from math import cos, sin
27 from numpy import dot, inner, sqrt, transpose
28
29
30 from lib.compat import norm
31 from lib.linear_algebra.kronecker_product import transpose_23
32 from lib.geometry.rotations import tilt_torsion_to_R
33
34
36 """Convert the frame order matrix (daeg) to the rotational superoperator.
37
38 @param daeg: The second degree frame order matrix, daeg. This must be in the Kronecker product layout.
39 @type daeg: numpy 9D, rank-2 array or numpy 3D, rank-4 array
40 @param Rsuper: The rotational superoperator structure to be populated.
41 @type Rsuper: numpy 5D, rank-2 array
42 """
43
44
45 transpose_23(daeg)
46
47
48 orig_shape = daeg.shape
49 daeg.shape = (3, 3, 3, 3)
50
51
52 Rsuper[0, 0] = daeg[0, 0, 0, 0] - daeg[2, 0, 2, 0]
53 Rsuper[1, 0] = daeg[0, 1, 0, 1] - daeg[2, 1, 2, 1]
54 Rsuper[2, 0] = daeg[0, 0, 0, 1] - daeg[2, 0, 2, 1]
55 Rsuper[3, 0] = daeg[0, 0, 0, 2] - daeg[2, 0, 2, 2]
56 Rsuper[4, 0] = daeg[0, 1, 0, 2] - daeg[2, 1, 2, 2]
57
58
59 Rsuper[0, 1] = daeg[1, 0, 1, 0] - daeg[2, 0, 2, 0]
60 Rsuper[1, 1] = daeg[1, 1, 1, 1] - daeg[2, 1, 2, 1]
61 Rsuper[2, 1] = daeg[1, 0, 1, 1] - daeg[2, 0, 2, 1]
62 Rsuper[3, 1] = daeg[1, 0, 1, 2] - daeg[2, 0, 2, 2]
63 Rsuper[4, 1] = daeg[1, 1, 1, 2] - daeg[2, 1, 2, 2]
64
65
66 Rsuper[0, 2] = daeg[0, 0, 1, 0] + daeg[1, 0, 0, 0]
67 Rsuper[1, 2] = daeg[0, 1, 1, 1] + daeg[1, 1, 0, 1]
68 Rsuper[2, 2] = daeg[0, 0, 1, 1] + daeg[1, 0, 0, 1]
69 Rsuper[3, 2] = daeg[0, 0, 1, 2] + daeg[1, 0, 0, 2]
70 Rsuper[4, 2] = daeg[0, 1, 1, 2] + daeg[1, 1, 0, 2]
71
72
73 Rsuper[0, 3] = daeg[0, 0, 2, 0] + daeg[2, 0, 0, 0]
74 Rsuper[1, 3] = daeg[0, 1, 2, 1] + daeg[2, 1, 0, 1]
75 Rsuper[2, 3] = daeg[0, 0, 2, 1] + daeg[2, 0, 0, 1]
76 Rsuper[3, 3] = daeg[0, 0, 2, 2] + daeg[2, 0, 0, 2]
77 Rsuper[4, 3] = daeg[0, 1, 2, 2] + daeg[2, 1, 0, 2]
78
79
80 Rsuper[0, 4] = daeg[1, 0, 2, 0] + daeg[2, 0, 1, 0]
81 Rsuper[1, 4] = daeg[1, 1, 2, 1] + daeg[2, 1, 1, 1]
82 Rsuper[2, 4] = daeg[1, 0, 2, 1] + daeg[2, 0, 1, 1]
83 Rsuper[3, 4] = daeg[1, 0, 2, 2] + daeg[2, 0, 1, 2]
84 Rsuper[4, 4] = daeg[1, 1, 2, 2] + daeg[2, 1, 1, 2]
85
86
87 daeg.shape = orig_shape
88
89
90 transpose_23(daeg)
91
92
93 -def pcs_pivot_motion_full(theta_i, phi_i, sigma_i, r_pivot_atom, r_ln_pivot, A, R_eigen, RT_eigen, Ri_prime):
94 """Calculate the PCS value after a pivoted motion for the isotropic cone model.
95
96 @param theta_i: The half cone opening angle (polar angle).
97 @type theta_i: float
98 @param phi_i: The cone azimuthal angle.
99 @type phi_i: float
100 @param sigma_i: The torsion angle for state i.
101 @type sigma_i: float
102 @param r_pivot_atom: The pivot point to atom vector.
103 @type r_pivot_atom: numpy rank-1, 3D array
104 @param r_ln_pivot: The lanthanide position to pivot point vector.
105 @type r_ln_pivot: numpy rank-1, 3D array
106 @param A: The full alignment tensor of the non-moving domain.
107 @type A: numpy rank-2, 3D array
108 @param R_eigen: The eigenframe rotation matrix.
109 @type R_eigen: numpy rank-2, 3D array
110 @param RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations).
111 @type RT_eigen: numpy rank-2, 3D array
112 @param Ri_prime: The empty rotation matrix for the in-frame isotropic cone motion for state i.
113 @type Ri_prime: numpy rank-2, 3D array
114 @return: The PCS value for the changed position.
115 @rtype: float
116 """
117
118
119 c_theta = cos(theta_i)
120 s_theta = sin(theta_i)
121 c_phi = cos(phi_i)
122 s_phi = sin(phi_i)
123 c_sigma_phi = cos(sigma_i - phi_i)
124 s_sigma_phi = sin(sigma_i - phi_i)
125 c_phi_c_theta = c_phi * c_theta
126 s_phi_c_theta = s_phi * c_theta
127 Ri_prime[0, 0] = c_phi_c_theta*c_sigma_phi - s_phi*s_sigma_phi
128 Ri_prime[0, 1] = -c_phi_c_theta*s_sigma_phi - s_phi*c_sigma_phi
129 Ri_prime[0, 2] = c_phi*s_theta
130 Ri_prime[1, 0] = s_phi_c_theta*c_sigma_phi + c_phi*s_sigma_phi
131 Ri_prime[1, 1] = -s_phi_c_theta*s_sigma_phi + c_phi*c_sigma_phi
132 Ri_prime[1, 2] = s_phi*s_theta
133 Ri_prime[2, 0] = -s_theta*c_sigma_phi
134 Ri_prime[2, 1] = s_theta*s_sigma_phi
135 Ri_prime[2, 2] = c_theta
136
137
138 R_i = dot(R_eigen, dot(Ri_prime, RT_eigen))
139
140
141 vect = dot(R_i, r_pivot_atom) + r_ln_pivot
142
143
144 length = norm(vect)
145
146
147 proj = dot(vect, dot(A, vect))
148
149
150 pcs = proj / length**5 * s_theta
151
152
153 return pcs
154
155
156 -def pcs_pivot_motion_full_qrint(theta_i=None, phi_i=None, sigma_i=None, full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None, r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None, pcs_theta=None, pcs_theta_err=None, missing_pcs=None, error_flag=False):
157 """Calculate the PCS value after a pivoted motion for the isotropic cone model.
158
159 @keyword theta_i: The half cone opening angle (polar angle).
160 @type theta_i: float
161 @keyword phi_i: The cone azimuthal angle.
162 @type phi_i: float
163 @keyword sigma_i: The torsion angle for state i.
164 @type sigma_i: float
165 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor.
166 @type full_in_ref_frame: numpy rank-1 array
167 @keyword r_pivot_atom: The pivot point to atom vector.
168 @type r_pivot_atom: numpy rank-2, 3D array
169 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector.
170 @type r_pivot_atom_rev: numpy rank-2, 3D array
171 @keyword r_ln_pivot: The lanthanide position to pivot point vector.
172 @type r_ln_pivot: numpy rank-2, 3D array
173 @keyword A: The full alignment tensor of the non-moving domain.
174 @type A: numpy rank-2, 3D array
175 @keyword R_eigen: The eigenframe rotation matrix.
176 @type R_eigen: numpy rank-2, 3D array
177 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations).
178 @type RT_eigen: numpy rank-2, 3D array
179 @keyword Ri_prime: The empty rotation matrix for the in-frame isotropic cone motion for state i.
180 @type Ri_prime: numpy rank-2, 3D array
181 @keyword pcs_theta: The storage structure for the back-calculated PCS values.
182 @type pcs_theta: numpy rank-2 array
183 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors.
184 @type pcs_theta_err: numpy rank-2 array
185 @keyword missing_pcs: A structure used to indicate which PCS values are missing.
186 @type missing_pcs: numpy rank-2 array
187 @keyword error_flag: A flag which if True will cause the PCS errors to be estimated and stored in pcs_theta_err.
188 @type error_flag: bool
189 """
190
191
192 tilt_torsion_to_R(phi_i, theta_i, sigma_i, Ri_prime)
193
194
195 R_i = dot(R_eigen, dot(Ri_prime, RT_eigen))
196
197
198 rot_vect_rev = transpose(dot(R_i, r_pivot_atom_rev) + r_ln_pivot)
199 rot_vect = transpose(dot(R_i, r_pivot_atom) + r_ln_pivot)
200
201
202 length_rev = 1.0 / norm(rot_vect_rev, axis=1)**5
203 length = 1.0 / norm(rot_vect, axis=1)**5
204
205
206 for j in range(len(r_pivot_atom[0])):
207
208 for i in range(len(pcs_theta)):
209
210 if missing_pcs[i, j]:
211 continue
212
213
214 if full_in_ref_frame[i]:
215 proj = dot(rot_vect[j], dot(A[i], rot_vect[j]))
216 length_i = length[j]
217 else:
218 proj = dot(rot_vect_rev[j], dot(A[i], rot_vect_rev[j]))
219 length_i = length_rev[j]
220
221
222 pcs_theta[i, j] += proj * length_i
223
224
225 if error_flag:
226 pcs_theta_err[i, j] += (proj * length_i)**2
227
228
230 """Calculate the PCS value after a pivoted motion for the isotropic cone model.
231
232 @param theta_i: The half cone opening angle (polar angle).
233 @type theta_i: float
234 @param phi_i: The cone azimuthal angle.
235 @type phi_i: float
236 @param r_pivot_atom: The pivot point to atom vector.
237 @type r_pivot_atom: numpy rank-1, 3D array
238 @param r_ln_pivot: The lanthanide position to pivot point vector.
239 @type r_ln_pivot: numpy rank-1, 3D array
240 @param A: The full alignment tensor of the non-moving domain.
241 @type A: numpy rank-2, 3D array
242 @param R_eigen: The eigenframe rotation matrix.
243 @type R_eigen: numpy rank-2, 3D array
244 @param RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations).
245 @type RT_eigen: numpy rank-2, 3D array
246 @param Ri_prime: The empty rotation matrix for the in-frame isotropic cone motion for state i.
247 @type Ri_prime: numpy rank-2, 3D array
248 @return: The PCS value for the changed position.
249 @rtype: float
250 """
251
252
253 c_theta = cos(theta_i)
254 s_theta = sin(theta_i)
255 c_phi = cos(phi_i)
256 s_phi = sin(phi_i)
257 c_phi_c_theta = c_phi * c_theta
258 s_phi_c_theta = s_phi * c_theta
259 Ri_prime[0, 0] = c_phi_c_theta*c_phi + s_phi**2
260 Ri_prime[0, 1] = c_phi_c_theta*s_phi - c_phi*s_phi
261 Ri_prime[0, 2] = c_phi*s_theta
262 Ri_prime[1, 0] = s_phi_c_theta*c_phi - c_phi*s_phi
263 Ri_prime[1, 1] = s_phi_c_theta*s_phi + c_phi**2
264 Ri_prime[1, 2] = s_phi*s_theta
265 Ri_prime[2, 0] = -s_theta*c_phi
266 Ri_prime[2, 1] = -s_theta*s_phi
267 Ri_prime[2, 2] = c_theta
268
269
270 R_i = dot(R_eigen, dot(Ri_prime, RT_eigen))
271
272
273 vect = dot(R_i, r_pivot_atom) + r_ln_pivot
274
275
276 length = norm(vect)
277
278
279 proj = dot(vect, dot(A, vect))
280
281
282 pcs = proj / length**5 * s_theta
283
284
285 return pcs
286
287
288 -def pcs_pivot_motion_torsionless_qrint(theta_i=None, phi_i=None, full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None, r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None, pcs_theta=None, pcs_theta_err=None, missing_pcs=None, error_flag=False):
289 """Calculate the PCS value after a pivoted motion for the isotropic cone model.
290
291 @keyword theta_i: The half cone opening angle (polar angle).
292 @type theta_i: float
293 @keyword phi_i: The cone azimuthal angle.
294 @type phi_i: float
295 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor.
296 @type full_in_ref_frame: numpy rank-1 array
297 @keyword r_pivot_atom: The pivot point to atom vector.
298 @type r_pivot_atom: numpy rank-2, 3D array
299 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector.
300 @type r_pivot_atom_rev: numpy rank-2, 3D array
301 @keyword r_ln_pivot: The lanthanide position to pivot point vector.
302 @type r_ln_pivot: numpy rank-2, 3D array
303 @keyword A: The full alignment tensor of the non-moving domain.
304 @type A: numpy rank-2, 3D array
305 @keyword R_eigen: The eigenframe rotation matrix.
306 @type R_eigen: numpy rank-2, 3D array
307 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations).
308 @type RT_eigen: numpy rank-2, 3D array
309 @keyword Ri_prime: The empty rotation matrix for the in-frame isotropic cone motion for state i.
310 @type Ri_prime: numpy rank-2, 3D array
311 @keyword pcs_theta: The storage structure for the back-calculated PCS values.
312 @type pcs_theta: numpy rank-2 array
313 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors.
314 @type pcs_theta_err: numpy rank-2 array
315 @keyword missing_pcs: A structure used to indicate which PCS values are missing.
316 @type missing_pcs: numpy rank-2 array
317 @keyword error_flag: A flag which if True will cause the PCS errors to be estimated and stored in pcs_theta_err.
318 @type error_flag: bool
319 """
320
321
322 c_theta = cos(theta_i)
323 s_theta = sin(theta_i)
324 c_phi = cos(phi_i)
325 s_phi = sin(phi_i)
326 c_phi_c_theta = c_phi * c_theta
327 s_phi_c_theta = s_phi * c_theta
328 Ri_prime[0, 0] = c_phi_c_theta*c_phi + s_phi**2
329 Ri_prime[0, 1] = c_phi_c_theta*s_phi - c_phi*s_phi
330 Ri_prime[0, 2] = c_phi*s_theta
331 Ri_prime[1, 0] = s_phi_c_theta*c_phi - c_phi*s_phi
332 Ri_prime[1, 1] = s_phi_c_theta*s_phi + c_phi**2
333 Ri_prime[1, 2] = s_phi*s_theta
334 Ri_prime[2, 0] = -s_theta*c_phi
335 Ri_prime[2, 1] = -s_theta*s_phi
336 Ri_prime[2, 2] = c_theta
337
338
339 R_i = dot(R_eigen, dot(Ri_prime, RT_eigen))
340
341
342 rot_vect_rev = transpose(dot(R_i, r_pivot_atom_rev) + r_ln_pivot)
343 rot_vect = transpose(dot(R_i, r_pivot_atom) + r_ln_pivot)
344
345
346 length_rev = 1.0 / norm(rot_vect_rev, axis=1)**5
347 length = 1.0 / norm(rot_vect, axis=1)**5
348
349
350 for j in range(len(r_pivot_atom[0])):
351
352 for i in range(len(pcs_theta)):
353
354 if missing_pcs[i, j]:
355 continue
356
357
358 if full_in_ref_frame[i]:
359 proj = dot(rot_vect[j], dot(A[i], rot_vect[j]))
360 length_i = length[j]
361 else:
362 proj = dot(rot_vect_rev[j], dot(A[i], rot_vect_rev[j]))
363 length_i = length_rev[j]
364
365
366 pcs_theta[i, j] += proj * length_i
367
368
369 if error_flag:
370 pcs_theta_err[i, j] += (proj * length_i)**2
371
372
374 """Calculate the reduction in the alignment tensor caused by the Frame Order matrix.
375
376 This is both the forward rotation notation and Kronecker product arrangement.
377
378 @param D: The Frame Order matrix, 2nd degree to be populated.
379 @type D: numpy 9D, rank-2 array
380 @param A: The full alignment tensor in {Axx, Ayy, Axy, Axz, Ayz} notation.
381 @type A: numpy 5D, rank-1 array
382 @param red_tensor: The structure in {Axx, Ayy, Axy, Axz, Ayz} notation to place the reduced
383 alignment tensor.
384 @type red_tensor: numpy 5D, rank-1 array
385 """
386
387
388 red_tensor[0] = (D[0, 0] - D[8, 0])*A[0]
389 red_tensor[0] = red_tensor[0] + (D[4, 0] - D[8, 0])*A[1]
390 red_tensor[0] = red_tensor[0] + (D[1, 0] + D[3, 0])*A[2]
391 red_tensor[0] = red_tensor[0] + (D[2, 0] + D[6, 0])*A[3]
392 red_tensor[0] = red_tensor[0] + (D[5, 0] + D[7, 0])*A[4]
393
394
395 red_tensor[1] = (D[0, 4] - D[8, 4])*A[0]
396 red_tensor[1] = red_tensor[1] + (D[4, 4] - D[8, 4])*A[1]
397 red_tensor[1] = red_tensor[1] + (D[1, 4] + D[3, 4])*A[2]
398 red_tensor[1] = red_tensor[1] + (D[2, 4] + D[6, 4])*A[3]
399 red_tensor[1] = red_tensor[1] + (D[5, 4] + D[7, 4])*A[4]
400
401
402 red_tensor[2] = (D[0, 1] - D[8, 1])*A[0]
403 red_tensor[2] = red_tensor[2] + (D[4, 1] - D[8, 1])*A[1]
404 red_tensor[2] = red_tensor[2] + (D[1, 1] + D[3, 1])*A[2]
405 red_tensor[2] = red_tensor[2] + (D[2, 1] + D[6, 1])*A[3]
406 red_tensor[2] = red_tensor[2] + (D[5, 1] + D[7, 1])*A[4]
407
408
409 red_tensor[3] = (D[0, 2] - D[8, 2])*A[0]
410 red_tensor[3] = red_tensor[3] + (D[4, 2] - D[8, 2])*A[1]
411 red_tensor[3] = red_tensor[3] + (D[1, 2] + D[3, 2])*A[2]
412 red_tensor[3] = red_tensor[3] + (D[2, 2] + D[6, 2])*A[3]
413 red_tensor[3] = red_tensor[3] + (D[5, 2] + D[7, 2])*A[4]
414
415
416 red_tensor[4] = (D[0, 5] - D[8, 5])*A[0]
417 red_tensor[4] = red_tensor[4] + (D[4, 5] - D[8, 5])*A[1]
418 red_tensor[4] = red_tensor[4] + (D[1, 5] + D[3, 5])*A[2]
419 red_tensor[4] = red_tensor[4] + (D[2, 5] + D[6, 5])*A[3]
420 red_tensor[4] = red_tensor[4] + (D[5, 5] + D[7, 5])*A[4]
421
422
424 """Calculate the reduction in the alignment tensor caused by the Frame Order matrix.
425
426 This is both the forward rotation notation and Kronecker product arrangement. This simplification is due to the symmetry in motion of the pseudo-elliptic and isotropic cones. All element of the frame order matrix where an index appears only once are zero.
427
428 @param D: The Frame Order matrix, 2nd degree to be populated.
429 @type D: numpy 9D, rank-2 array
430 @param A: The full alignment tensor in {Axx, Ayy, Axy, Axz, Ayz} notation.
431 @type A: numpy 5D, rank-1 array
432 @param red_tensor: The structure in {Axx, Ayy, Axy, Axz, Ayz} notation to place the reduced
433 alignment tensor.
434 @type red_tensor: numpy 5D, rank-1 array
435 """
436
437
438 red_tensor[0] = (D[0, 0] - D[8, 0])*A[0] + (D[4, 0] - D[8, 0])*A[1]
439 red_tensor[1] = (D[0, 4] - D[8, 4])*A[0] + (D[4, 4] - D[8, 4])*A[1]
440 red_tensor[2] = (D[1, 1] + D[3, 1])*A[2]
441 red_tensor[3] = (D[2, 2] + D[6, 2])*A[3]
442 red_tensor[4] = (D[5, 5] + D[7, 5])*A[4]
443
444
446 """Rotate the given frame order matrix.
447
448 It is assumed that the frame order matrix is in the Kronecker product form.
449
450
451 @param matrix: The Frame Order matrix, 2nd degree to be populated.
452 @type matrix: numpy 9D, rank-2 array
453 @param Rx2_eigen: The Kronecker product of the eigenframe rotation matrix with itself.
454 @type Rx2_eigen: numpy 9D, rank-2 array
455 """
456
457
458 matrix_rot = dot(Rx2_eigen, dot(matrix, transpose(Rx2_eigen)))
459
460
461 return matrix_rot
462
463
465 """A data container stored in the memo objects for use by the Result_command class."""
466