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, pi, sin, sqrt
27 try:
28 from scipy.integrate import dblquad, quad
29 except ImportError:
30 pass
31
32
33 from lib.geometry.pec import pec
34 from lib.frame_order.matrix_ops import pcs_pivot_motion_torsionless, pcs_pivot_motion_torsionless_qrint, rotate_daeg
35 from lib.frame_order.pseudo_ellipse import tmax_pseudo_ellipse
36
37
39 """Generate the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse.
40
41 @param matrix: The Frame Order matrix, 2nd degree to be populated.
42 @type matrix: numpy 9D, rank-2 array
43 @param Rx2_eigen: The Kronecker product of the eigenframe rotation matrix with itself.
44 @type Rx2_eigen: numpy 9D, rank-2 array
45 @param theta_x: The cone opening angle along x.
46 @type theta_x: float
47 @param theta_y: The cone opening angle along y.
48 @type theta_y: float
49 """
50
51
52 fact = 1.0 / (6.0 * pec(theta_x, theta_y))
53
54
55 matrix[0, 0] = fact * (6.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_00, -pi, pi, args=(theta_x, theta_y), full_output=1)[0])
56 matrix[1, 1] = fact * (2.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_11, -pi, pi, args=(theta_x, theta_y), full_output=1)[0])
57 matrix[2, 2] = fact * (5.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_22, -pi, pi, args=(theta_x, theta_y), full_output=1)[0])
58 matrix[3, 3] = matrix[1, 1]
59 matrix[4, 4] = fact * (6.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_44, -pi, pi, args=(theta_x, theta_y), full_output=1)[0])
60 matrix[5, 5] = fact * (5.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_55, -pi, pi, args=(theta_x, theta_y), full_output=1)[0])
61 matrix[6, 6] = matrix[2, 2]
62 matrix[7, 7] = matrix[5, 5]
63 matrix[8, 8] = fact * quad(part_int_daeg2_pseudo_ellipse_torsionless_88, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]
64
65
66 matrix[0, 4] = matrix[4, 0] = fact * (2.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_04, -pi, pi, args=(theta_x, theta_y), full_output=1)[0])
67 matrix[0, 8] = matrix[8, 0] = fact * (4.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_08, -pi, pi, args=(theta_x, theta_y), full_output=1)[0])
68 matrix[4, 8] = matrix[8, 4] = fact * (4.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_48, -pi, pi, args=(theta_x, theta_y), full_output=1)[0])
69
70
71 matrix[1, 3] = matrix[3, 1] = matrix[0, 4]
72 matrix[2, 6] = matrix[6, 2] = -matrix[0, 8]
73 matrix[5, 7] = matrix[7, 5] = -matrix[4, 8]
74
75
76 return rotate_daeg(matrix, Rx2_eigen)
77
78
80 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse.
81
82 @param phi: The azimuthal tilt-torsion angle.
83 @type phi: float
84 @param x: The cone opening angle along x.
85 @type x: float
86 @param y: The cone opening angle along y.
87 @type y: float
88 @return: The theta partial integral.
89 @rtype: float
90 """
91
92
93 tmax = tmax_pseudo_ellipse(phi, x, y)
94
95
96 return (2*cos(phi)**4*cos(tmax) + 6*cos(phi)**2*sin(phi)**2)*sin(tmax)**2 - (6*sin(phi)**4 + 2*cos(phi)**4)*cos(tmax)
97
98
100 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse.
101
102 @param phi: The azimuthal tilt-torsion angle.
103 @type phi: float
104 @param x: The cone opening angle along x.
105 @type x: float
106 @param y: The cone opening angle along y.
107 @type y: float
108 @return: The theta partial integral.
109 @rtype: float
110 """
111
112
113 tmax = tmax_pseudo_ellipse(phi, x, y)
114
115
116 return (2*cos(phi)**2*sin(phi)**2*cos(tmax) - 6*cos(phi)**2*sin(phi)**2)*sin(tmax)**2 - 8*cos(phi)**2*sin(phi)**2*cos(tmax)
117
118
120 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse.
121
122 @param phi: The azimuthal tilt-torsion angle.
123 @type phi: float
124 @param x: The cone opening angle along x.
125 @type x: float
126 @param y: The cone opening angle along y.
127 @type y: float
128 @return: The theta partial integral.
129 @rtype: float
130 """
131
132
133 tmax = tmax_pseudo_ellipse(phi, x, y)
134
135
136 return 2*cos(phi)**2*cos(tmax)**3 - 6*cos(phi)**2*cos(tmax)
137
138
140 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse.
141
142 @param phi: The azimuthal tilt-torsion angle.
143 @type phi: float
144 @param x: The cone opening angle along x.
145 @type x: float
146 @param y: The cone opening angle along y.
147 @type y: float
148 @return: The theta partial integral.
149 @rtype: float
150 """
151
152
153 tmax = tmax_pseudo_ellipse(phi, x, y)
154
155
156 return (2*cos(phi)**2*sin(phi)**2*cos(tmax) + 3*sin(phi)**4 + 3*cos(phi)**4)*sin(tmax)**2 - 8*cos(phi)**2*sin(phi)**2*cos(tmax)
157
158
160 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse.
161
162 @param phi: The azimuthal tilt-torsion angle.
163 @type phi: float
164 @param x: The cone opening angle along x.
165 @type x: float
166 @param y: The cone opening angle along y.
167 @type y: float
168 @return: The theta partial integral.
169 @rtype: float
170 """
171
172
173 tmax = tmax_pseudo_ellipse(phi, x, y)
174
175
176 return (2*sin(phi)**2 - 2)*cos(tmax)**3 - 3*sin(phi)**2*cos(tmax)**2
177
178
180 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse.
181
182 @param phi: The azimuthal tilt-torsion angle.
183 @type phi: float
184 @param x: The cone opening angle along x.
185 @type x: float
186 @param y: The cone opening angle along y.
187 @type y: float
188 @return: The theta partial integral.
189 @rtype: float
190 """
191
192
193 tmax = tmax_pseudo_ellipse(phi, x, y)
194
195
196 return (2*sin(phi)**4*cos(tmax) + 6*cos(phi)**2*sin(phi)**2)*sin(tmax)**2 - (2*sin(phi)**4 + 6*cos(phi)**4)*cos(tmax)
197
198
200 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse.
201
202 @param phi: The azimuthal tilt-torsion angle.
203 @type phi: float
204 @param x: The cone opening angle along x.
205 @type x: float
206 @param y: The cone opening angle along y.
207 @type y: float
208 @return: The theta partial integral.
209 @rtype: float
210 """
211
212
213 tmax = tmax_pseudo_ellipse(phi, x, y)
214
215
216 return 2*sin(phi)**2*cos(tmax)**3 - 6*sin(phi)**2*cos(tmax)
217
218
220 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse.
221
222 @param phi: The azimuthal tilt-torsion angle.
223 @type phi: float
224 @param x: The cone opening angle along x.
225 @type x: float
226 @param y: The cone opening angle along y.
227 @type y: float
228 @return: The theta partial integral.
229 @rtype: float
230 """
231
232
233 tmax = tmax_pseudo_ellipse(phi, x, y)
234
235
236 return (2*cos(phi)**2 - 2)*cos(tmax)**3 - 3*cos(phi)**2*cos(tmax)**2
237
238
240 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse.
241
242 @param phi: The azimuthal tilt-torsion angle.
243 @type phi: float
244 @param x: The cone opening angle along x.
245 @type x: float
246 @param y: The cone opening angle along y.
247 @type y: float
248 @return: The theta partial integral.
249 @rtype: float
250 """
251
252
253 tmax = tmax_pseudo_ellipse(phi, x, y)
254
255
256 return 2 - 2*cos(tmax)**3
257
258
260 """Determine the averaged PCS value via numerical integration.
261
262 @keyword theta_x: The x-axis half cone angle.
263 @type theta_x: float
264 @keyword theta_y: The y-axis half cone angle.
265 @type theta_y: float
266 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units).
267 @type c: float
268 @keyword r_pivot_atom: The pivot point to atom vector.
269 @type r_pivot_atom: numpy rank-1, 3D array
270 @keyword r_ln_pivot: The lanthanide position to pivot point vector.
271 @type r_ln_pivot: numpy rank-1, 3D array
272 @keyword A: The full alignment tensor of the non-moving domain.
273 @type A: numpy rank-2, 3D array
274 @keyword R_eigen: The eigenframe rotation matrix.
275 @type R_eigen: numpy rank-2, 3D array
276 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations).
277 @type RT_eigen: numpy rank-2, 3D array
278 @keyword Ri_prime: The empty rotation matrix for the in-frame isotropic cone motion, used to calculate the PCS for each state i in the numerical integration.
279 @type Ri_prime: numpy rank-2, 3D array
280 @return: The averaged PCS value.
281 @rtype: float
282 """
283
284 def pseudo_ellipse(phi):
285 """The pseudo-ellipse wrapper formula."""
286
287 return tmax_pseudo_ellipse(phi, theta_x, theta_y)
288
289
290 result = dblquad(pcs_pivot_motion_torsionless, -pi, pi, lambda phi: 0.0, pseudo_ellipse, args=(r_pivot_atom, r_ln_pivot, A, R_eigen, RT_eigen, Ri_prime))
291
292
293 SA = pec(theta_x, theta_y)
294
295
296 return c * result[0] / SA
297
298
299 -def pcs_numeric_int_pseudo_ellipse_torsionless_qrint(points=None, theta_x=None, theta_y=None, c=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):
300 """Determine the averaged PCS value via numerical integration.
301
302 @keyword points: The Sobol points in the torsion-tilt angle space.
303 @type points: numpy rank-2, 3D array
304 @keyword theta_x: The x-axis half cone angle.
305 @type theta_x: float
306 @keyword theta_y: The y-axis half cone angle.
307 @type theta_y: float
308 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units).
309 @type c: numpy rank-1 array
310 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor.
311 @type full_in_ref_frame: numpy rank-1 array
312 @keyword r_pivot_atom: The pivot point to atom vector.
313 @type r_pivot_atom: numpy rank-2, 3D array
314 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector.
315 @type r_pivot_atom_rev: numpy rank-2, 3D array
316 @keyword r_ln_pivot: The lanthanide position to pivot point vector.
317 @type r_ln_pivot: numpy rank-2, 3D array
318 @keyword A: The full alignment tensor of the non-moving domain.
319 @type A: numpy rank-2, 3D array
320 @keyword R_eigen: The eigenframe rotation matrix.
321 @type R_eigen: numpy rank-2, 3D array
322 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations).
323 @type RT_eigen: numpy rank-2, 3D array
324 @keyword Ri_prime: The empty rotation matrix for the in-frame isotropic cone motion, used to calculate the PCS for each state i in the numerical integration.
325 @type Ri_prime: numpy rank-2, 3D array
326 @keyword pcs_theta: The storage structure for the back-calculated PCS values.
327 @type pcs_theta: numpy rank-2 array
328 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors.
329 @type pcs_theta_err: numpy rank-2 array
330 @keyword missing_pcs: A structure used to indicate which PCS values are missing.
331 @type missing_pcs: numpy rank-2 array
332 @keyword error_flag: A flag which if True will cause the PCS errors to be estimated and stored in pcs_theta_err.
333 @type error_flag: bool
334 """
335
336
337 for i in range(len(pcs_theta)):
338 for j in range(len(pcs_theta[i])):
339 pcs_theta[i, j] = 0.0
340 pcs_theta_err[i, j] = 0.0
341
342
343 num = 0
344 for i in range(len(points)):
345
346 theta, phi = points[i]
347
348
349 theta_max = tmax_pseudo_ellipse(phi, theta_x, theta_y)
350
351
352 if theta > theta_max:
353 continue
354
355
356 pcs_pivot_motion_torsionless_qrint(theta_i=theta, phi_i=phi, full_in_ref_frame=full_in_ref_frame, r_pivot_atom=r_pivot_atom, r_pivot_atom_rev=r_pivot_atom_rev, r_ln_pivot=r_ln_pivot, A=A, R_eigen=R_eigen, RT_eigen=RT_eigen, Ri_prime=Ri_prime, pcs_theta=pcs_theta, pcs_theta_err=pcs_theta_err, missing_pcs=missing_pcs)
357
358
359 num += 1
360
361
362 for i in range(len(pcs_theta)):
363 for j in range(len(pcs_theta[i])):
364
365 pcs_theta[i, j] = c[i] * pcs_theta[i, j] / float(num)
366
367
368 if error_flag:
369 pcs_theta_err[i, j] = abs(pcs_theta_err[i, j] / float(num) - pcs_theta[i, j]**2) / float(num)
370 pcs_theta_err[i, j] = c[i] * sqrt(pcs_theta_err[i, j])
371 print("%8.3f +/- %-8.3f" % (pcs_theta[i, j]*1e6, pcs_theta_err[i, j]*1e6))
372