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, transpose
28 from numpy.linalg import norm
29
30
31 from lib.compat import norm
32 from lib.linear_algebra.kronecker_product import transpose_23
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_qr_int(full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None, r_ln_pivot=None, A=None, Ri=None, pcs_theta=None, pcs_theta_err=None, missing_pcs=None):
94 """Calculate the PCS value after a pivoted motion for the isotropic cone model.
95
96 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor.
97 @type full_in_ref_frame: numpy rank-1 array
98 @keyword r_pivot_atom: The pivot point to atom vector.
99 @type r_pivot_atom: numpy rank-2, 3D array
100 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector.
101 @type r_pivot_atom_rev: numpy rank-2, 3D array
102 @keyword r_ln_pivot: The lanthanide position to pivot point vector.
103 @type r_ln_pivot: numpy rank-2, 3D array
104 @keyword A: The full alignment tensor of the non-moving domain.
105 @type A: numpy rank-2, 3D array
106 @keyword Ri: The frame-shifted, pre-calculated rotation matrix for state i.
107 @type Ri: numpy rank-2, 3D array
108 @keyword pcs_theta: The storage structure for the back-calculated PCS values.
109 @type pcs_theta: numpy rank-2 array
110 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors.
111 @type pcs_theta_err: numpy rank-2 array
112 @keyword missing_pcs: A structure used to indicate which PCS values are missing.
113 @type missing_pcs: numpy rank-2 array
114 """
115
116
117 rot_vect = dot(r_pivot_atom, Ri) + r_ln_pivot
118
119
120 length = 1.0 / norm(rot_vect, axis=1)**5
121
122
123 if min(full_in_ref_frame) == 0:
124 rot_vect_rev = dot(r_pivot_atom_rev, Ri) + r_ln_pivot
125 length_rev = 1.0 / norm(rot_vect_rev, axis=1)**5
126
127
128 for j in range(len(r_pivot_atom[:, 0])):
129
130 for i in range(len(pcs_theta)):
131
132 if missing_pcs[i, j]:
133 continue
134
135
136 if full_in_ref_frame[i]:
137 proj = dot(rot_vect[j], dot(A[i], rot_vect[j]))
138 length_i = length[j]
139 else:
140 proj = dot(rot_vect_rev[j], dot(A[i], rot_vect_rev[j]))
141 length_i = length_rev[j]
142
143
144 pcs_theta[i, j] += proj * length_i
145
146
148 """Calculate the PCS value after a pivoted motion for the isotropic cone model.
149
150 @param theta_i: The half cone opening angle (polar angle).
151 @type theta_i: float
152 @param phi_i: The cone azimuthal angle.
153 @type phi_i: float
154 @param sigma_i: The torsion angle for state i.
155 @type sigma_i: float
156 @param r_pivot_atom: The pivot point to atom vector.
157 @type r_pivot_atom: numpy rank-1, 3D array
158 @param r_ln_pivot: The lanthanide position to pivot point vector.
159 @type r_ln_pivot: numpy rank-1, 3D array
160 @param A: The full alignment tensor of the non-moving domain.
161 @type A: numpy rank-2, 3D array
162 @param R_eigen: The eigenframe rotation matrix.
163 @type R_eigen: numpy rank-2, 3D array
164 @param RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations).
165 @type RT_eigen: numpy rank-2, 3D array
166 @param Ri_prime: The empty rotation matrix for the in-frame isotropic cone motion for state i.
167 @type Ri_prime: numpy rank-2, 3D array
168 @return: The PCS value for the changed position.
169 @rtype: float
170 """
171
172
173 c_theta = cos(theta_i)
174 s_theta = sin(theta_i)
175 c_phi = cos(phi_i)
176 s_phi = sin(phi_i)
177 c_sigma_phi = cos(sigma_i - phi_i)
178 s_sigma_phi = sin(sigma_i - phi_i)
179 c_phi_c_theta = c_phi * c_theta
180 s_phi_c_theta = s_phi * c_theta
181 Ri_prime[0, 0] = c_phi_c_theta*c_sigma_phi - s_phi*s_sigma_phi
182 Ri_prime[0, 1] = -c_phi_c_theta*s_sigma_phi - s_phi*c_sigma_phi
183 Ri_prime[0, 2] = c_phi*s_theta
184 Ri_prime[1, 0] = s_phi_c_theta*c_sigma_phi + c_phi*s_sigma_phi
185 Ri_prime[1, 1] = -s_phi_c_theta*s_sigma_phi + c_phi*c_sigma_phi
186 Ri_prime[1, 2] = s_phi*s_theta
187 Ri_prime[2, 0] = -s_theta*c_sigma_phi
188 Ri_prime[2, 1] = s_theta*s_sigma_phi
189 Ri_prime[2, 2] = c_theta
190
191
192 R_i = dot(R_eigen, dot(Ri_prime, RT_eigen))
193
194
195 vect = dot(R_i, r_pivot_atom) + r_ln_pivot
196
197
198 length = norm(vect)
199
200
201 proj = dot(vect, dot(A, vect))
202
203
204 pcs = proj / length**5 * s_theta
205
206
207 return pcs
208
209
210 -def pcs_pivot_motion_torsionless_qr_int(full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None, r_ln_pivot=None, A=None, Ri=None, pcs_theta=None, pcs_theta_err=None, missing_pcs=None):
211 """Calculate the PCS value after a pivoted motion for the isotropic cone model.
212
213 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor.
214 @type full_in_ref_frame: numpy rank-1 array
215 @keyword r_pivot_atom: The pivot point to atom vector.
216 @type r_pivot_atom: numpy rank-2, 3D array
217 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector.
218 @type r_pivot_atom_rev: numpy rank-2, 3D array
219 @keyword r_ln_pivot: The lanthanide position to pivot point vector.
220 @type r_ln_pivot: numpy rank-2, 3D array
221 @keyword A: The full alignment tensor of the non-moving domain.
222 @type A: numpy rank-2, 3D array
223 @keyword Ri: The frame-shifted, pre-calculated rotation matrix for state i.
224 @type Ri: numpy rank-2, 3D array
225 @keyword pcs_theta: The storage structure for the back-calculated PCS values.
226 @type pcs_theta: numpy rank-2 array
227 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors.
228 @type pcs_theta_err: numpy rank-2 array
229 @keyword missing_pcs: A structure used to indicate which PCS values are missing.
230 @type missing_pcs: numpy rank-2 array
231 """
232
233
234 rot_vect = dot(r_pivot_atom, Ri) + r_ln_pivot
235
236
237 length = 1.0 / norm(rot_vect, axis=1)**5
238
239
240 if min(full_in_ref_frame) == 0:
241 rot_vect_rev = dot(r_pivot_atom_rev, Ri) + r_ln_pivot
242 length_rev = 1.0 / norm(rot_vect_rev, axis=1)**5
243
244
245 for j in range(len(r_pivot_atom[:, 0])):
246
247 for i in range(len(pcs_theta)):
248
249 if missing_pcs[i, j]:
250 continue
251
252
253 if full_in_ref_frame[i]:
254 proj = dot(rot_vect[j], dot(A[i], rot_vect[j]))
255 length_i = length[j]
256 else:
257 proj = dot(rot_vect_rev[j], dot(A[i], rot_vect_rev[j]))
258 length_i = length_rev[j]
259
260
261 pcs_theta[i, j] += proj * length_i
262
263
265 """Calculate the PCS value after a pivoted motion for the isotropic cone model.
266
267 @param theta_i: The half cone opening angle (polar angle).
268 @type theta_i: float
269 @param phi_i: The cone azimuthal angle.
270 @type phi_i: float
271 @param r_pivot_atom: The pivot point to atom vector.
272 @type r_pivot_atom: numpy rank-1, 3D array
273 @param r_ln_pivot: The lanthanide position to pivot point vector.
274 @type r_ln_pivot: numpy rank-1, 3D array
275 @param A: The full alignment tensor of the non-moving domain.
276 @type A: numpy rank-2, 3D array
277 @param R_eigen: The eigenframe rotation matrix.
278 @type R_eigen: numpy rank-2, 3D array
279 @param RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations).
280 @type RT_eigen: numpy rank-2, 3D array
281 @param Ri_prime: The empty rotation matrix for the in-frame isotropic cone motion for state i.
282 @type Ri_prime: numpy rank-2, 3D array
283 @return: The PCS value for the changed position.
284 @rtype: float
285 """
286
287
288 c_theta = cos(theta_i)
289 s_theta = sin(theta_i)
290 c_phi = cos(phi_i)
291 s_phi = sin(phi_i)
292 c_phi_c_theta = c_phi * c_theta
293 s_phi_c_theta = s_phi * c_theta
294 Ri_prime[0, 0] = c_phi_c_theta*c_phi + s_phi**2
295 Ri_prime[0, 1] = c_phi_c_theta*s_phi - c_phi*s_phi
296 Ri_prime[0, 2] = c_phi*s_theta
297 Ri_prime[1, 0] = s_phi_c_theta*c_phi - c_phi*s_phi
298 Ri_prime[1, 1] = s_phi_c_theta*s_phi + c_phi**2
299 Ri_prime[1, 2] = s_phi*s_theta
300 Ri_prime[2, 0] = -s_theta*c_phi
301 Ri_prime[2, 1] = -s_theta*s_phi
302 Ri_prime[2, 2] = c_theta
303
304
305 R_i = dot(R_eigen, dot(Ri_prime, RT_eigen))
306
307
308 vect = dot(R_i, r_pivot_atom) + r_ln_pivot
309
310
311 length = norm(vect)
312
313
314 proj = dot(vect, dot(A, vect))
315
316
317 pcs = proj / length**5 * s_theta
318
319
320 return pcs
321
322
324 """Calculate the reduction in the alignment tensor caused by the Frame Order matrix.
325
326 This is both the forward rotation notation and Kronecker product arrangement.
327
328 @param D: The Frame Order matrix, 2nd degree to be populated.
329 @type D: numpy 9D, rank-2 array
330 @param A: The full alignment tensor in {Axx, Ayy, Axy, Axz, Ayz} notation.
331 @type A: numpy 5D, rank-1 array
332 @param red_tensor: The structure in {Axx, Ayy, Axy, Axz, Ayz} notation to place the reduced
333 alignment tensor.
334 @type red_tensor: numpy 5D, rank-1 array
335 """
336
337
338 red_tensor[0] = (D[0, 0] - D[8, 0])*A[0]
339 red_tensor[0] = red_tensor[0] + (D[4, 0] - D[8, 0])*A[1]
340 red_tensor[0] = red_tensor[0] + (D[1, 0] + D[3, 0])*A[2]
341 red_tensor[0] = red_tensor[0] + (D[2, 0] + D[6, 0])*A[3]
342 red_tensor[0] = red_tensor[0] + (D[5, 0] + D[7, 0])*A[4]
343
344
345 red_tensor[1] = (D[0, 4] - D[8, 4])*A[0]
346 red_tensor[1] = red_tensor[1] + (D[4, 4] - D[8, 4])*A[1]
347 red_tensor[1] = red_tensor[1] + (D[1, 4] + D[3, 4])*A[2]
348 red_tensor[1] = red_tensor[1] + (D[2, 4] + D[6, 4])*A[3]
349 red_tensor[1] = red_tensor[1] + (D[5, 4] + D[7, 4])*A[4]
350
351
352 red_tensor[2] = (D[0, 1] - D[8, 1])*A[0]
353 red_tensor[2] = red_tensor[2] + (D[4, 1] - D[8, 1])*A[1]
354 red_tensor[2] = red_tensor[2] + (D[1, 1] + D[3, 1])*A[2]
355 red_tensor[2] = red_tensor[2] + (D[2, 1] + D[6, 1])*A[3]
356 red_tensor[2] = red_tensor[2] + (D[5, 1] + D[7, 1])*A[4]
357
358
359 red_tensor[3] = (D[0, 2] - D[8, 2])*A[0]
360 red_tensor[3] = red_tensor[3] + (D[4, 2] - D[8, 2])*A[1]
361 red_tensor[3] = red_tensor[3] + (D[1, 2] + D[3, 2])*A[2]
362 red_tensor[3] = red_tensor[3] + (D[2, 2] + D[6, 2])*A[3]
363 red_tensor[3] = red_tensor[3] + (D[5, 2] + D[7, 2])*A[4]
364
365
366 red_tensor[4] = (D[0, 5] - D[8, 5])*A[0]
367 red_tensor[4] = red_tensor[4] + (D[4, 5] - D[8, 5])*A[1]
368 red_tensor[4] = red_tensor[4] + (D[1, 5] + D[3, 5])*A[2]
369 red_tensor[4] = red_tensor[4] + (D[2, 5] + D[6, 5])*A[3]
370 red_tensor[4] = red_tensor[4] + (D[5, 5] + D[7, 5])*A[4]
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. 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.
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] + (D[4, 0] - D[8, 0])*A[1]
389 red_tensor[1] = (D[0, 4] - D[8, 4])*A[0] + (D[4, 4] - D[8, 4])*A[1]
390 red_tensor[2] = (D[1, 1] + D[3, 1])*A[2]
391 red_tensor[3] = (D[2, 2] + D[6, 2])*A[3]
392 red_tensor[4] = (D[5, 5] + D[7, 5])*A[4]
393
394
396 """Rotate the given frame order matrix.
397
398 It is assumed that the frame order matrix is in the Kronecker product form.
399
400
401 @param matrix: The Frame Order matrix to be populated. All degrees are handled.
402 @type matrix: numpy rank-2 array
403 @param R_eigen: The Kronecker product of the eigenframe rotation matrix with itself.
404 @type R_eigen: numpy rank-2 array
405 """
406
407
408 matrix_rot = dot(R_eigen, dot(matrix, transpose(R_eigen)))
409
410
411 return matrix_rot
412
413
415 """A data container stored in the memo objects for use by the Result_command class."""
416