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 import dep_check
27
28
29 from math import cos, pi, sin, sqrt
30 from numpy import dot, inner, sinc, transpose
31 from numpy.linalg import norm
32 if dep_check.scipy_module:
33 from scipy.integrate import quad
34
35
36 from lib.frame_order.matrix_ops import rotate_daeg
37
38
40 """Generate the rotated 2nd degree Frame Order matrix for the rotor model.
41
42 The cone axis is assumed to be parallel to the z-axis in the eigenframe.
43
44
45 @param matrix: The Frame Order matrix, 2nd degree to be populated.
46 @type matrix: numpy 9D, rank-2 array
47 @param Rx2_eigen: The Kronecker product of the eigenframe rotation matrix with itself.
48 @type Rx2_eigen: numpy 9D, rank-2 array
49 @param smax: The maximum torsion angle.
50 @type smax: float
51 """
52
53
54 for i in range(9):
55 for j in range(9):
56 matrix[i, j] = 0.0
57
58
59 sinc_smax = sinc(smax/pi)
60 sinc_2smax = sinc(2.0*smax/pi)
61
62
63 matrix[0, 0] = (sinc_2smax + 1.0) / 2.0
64 matrix[1, 1] = matrix[0, 0]
65 matrix[2, 2] = sinc_smax
66 matrix[3, 3] = matrix[0, 0]
67 matrix[4, 4] = matrix[0, 0]
68 matrix[5, 5] = matrix[2, 2]
69 matrix[6, 6] = matrix[2, 2]
70 matrix[7, 7] = matrix[2, 2]
71 matrix[8, 8] = 1.0
72
73
74 matrix[0, 4] = matrix[4, 0] = -(sinc_2smax - 1.0) / 2.0
75
76
77 matrix[1, 3] = matrix[3, 1] = -matrix[0, 4]
78
79
80 return rotate_daeg(matrix, Rx2_eigen)
81
82
83 -def pcs_numeric_int_rotor(sigma_max=None, c=None, r_pivot_atom=None, r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None):
84 """Determine the averaged PCS value via numerical integration.
85
86 @keyword sigma_max: The maximum rotor angle.
87 @type sigma_max: float
88 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units).
89 @type c: float
90 @keyword r_pivot_atom: The pivot point to atom vector.
91 @type r_pivot_atom: numpy rank-1, 3D array
92 @keyword r_ln_pivot: The lanthanide position to pivot point vector.
93 @type r_ln_pivot: numpy rank-1, 3D array
94 @keyword A: The full alignment tensor of the non-moving domain.
95 @type A: numpy rank-2, 3D array
96 @keyword R_eigen: The eigenframe rotation matrix.
97 @type R_eigen: numpy rank-2, 3D array
98 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations).
99 @type RT_eigen: numpy rank-2, 3D array
100 @keyword Ri_prime: The empty rotation matrix for the in-frame rotor motion, used to calculate the PCS for each state i in the numerical integration.
101 @type Ri_prime: numpy rank-2, 3D array
102 @return: The averaged PCS value.
103 @rtype: float
104 """
105
106
107 Ri_prime[0, 2] = 0.0
108 Ri_prime[1, 2] = 0.0
109 Ri_prime[2, 0] = 0.0
110 Ri_prime[2, 1] = 0.0
111 Ri_prime[2, 2] = 1.0
112
113
114 result = quad(pcs_pivot_motion_rotor, -sigma_max, sigma_max, args=(r_pivot_atom, r_ln_pivot, A, R_eigen, RT_eigen, Ri_prime))
115
116
117 SA = 2.0 * sigma_max
118
119
120 return c * result[0] / SA
121
122
123 -def pcs_numeric_int_rotor_qrint(points=None, sigma_max=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):
124 """Determine the averaged PCS value via numerical integration.
125
126 @keyword points: The Sobol points in the torsion-tilt angle space.
127 @type points: numpy rank-2, 3D array
128 @keyword sigma_max: The maximum rotor angle.
129 @type sigma_max: float
130 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units).
131 @type c: numpy rank-1 array
132 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor.
133 @type full_in_ref_frame: numpy rank-1 array
134 @keyword r_pivot_atom: The pivot point to atom vector.
135 @type r_pivot_atom: numpy rank-2, 3D array
136 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector.
137 @type r_pivot_atom_rev: numpy rank-2, 3D array
138 @keyword r_ln_pivot: The lanthanide position to pivot point vector.
139 @type r_ln_pivot: numpy rank-2, 3D array
140 @keyword A: The full alignment tensor of the non-moving domain.
141 @type A: numpy rank-2, 3D array
142 @keyword R_eigen: The eigenframe rotation matrix.
143 @type R_eigen: numpy rank-2, 3D array
144 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations).
145 @type RT_eigen: numpy rank-2, 3D array
146 @keyword Ri_prime: The empty rotation matrix for the in-frame rotor motion, used to calculate the PCS for each state i in the numerical integration.
147 @type Ri_prime: numpy rank-2, 3D array
148 @keyword pcs_theta: The storage structure for the back-calculated PCS values.
149 @type pcs_theta: numpy rank-2 array
150 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors.
151 @type pcs_theta_err: numpy rank-2 array
152 @keyword missing_pcs: A structure used to indicate which PCS values are missing.
153 @type missing_pcs: numpy rank-2 array
154 @keyword error_flag: A flag which if True will cause the PCS errors to be estimated and stored in pcs_theta_err.
155 @type error_flag: bool
156 """
157
158
159 for i in range(len(pcs_theta)):
160 for j in range(len(pcs_theta[i])):
161 pcs_theta[i, j] = 0.0
162 pcs_theta_err[i, j] = 0.0
163
164
165 num = 0
166 for i in range(len(points)):
167
168 sigma = points[i]
169
170
171 if sigma > sigma_max or sigma < -sigma_max:
172 continue
173
174
175 pcs_pivot_motion_rotor_qrint(sigma_i=sigma, 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)
176
177
178 num += 1
179
180
181 for i in range(len(pcs_theta)):
182 for j in range(len(pcs_theta[i])):
183
184 pcs_theta[i, j] = c[i] * pcs_theta[i, j] / float(num)
185
186
187 if error_flag:
188 pcs_theta_err[i, j] = abs(pcs_theta_err[i, j] / float(num) - pcs_theta[i, j]**2) / float(num)
189 pcs_theta_err[i, j] = c[i] * sqrt(pcs_theta_err[i, j])
190 print("%8.3f +/- %-8.3f" % (pcs_theta[i, j]*1e6, pcs_theta_err[i, j]*1e6))
191
192
194 """Calculate the PCS value after a pivoted motion for the rotor model.
195
196 @param sigma_i: The rotor angle for state i.
197 @type sigma_i: float
198 @param r_pivot_atom: The pivot point to atom vector.
199 @type r_pivot_atom: numpy rank-1, 3D array
200 @param r_ln_pivot: The lanthanide position to pivot point vector.
201 @type r_ln_pivot: numpy rank-1, 3D array
202 @param A: The full alignment tensor of the non-moving domain.
203 @type A: numpy rank-2, 3D array
204 @param R_eigen: The eigenframe rotation matrix.
205 @type R_eigen: numpy rank-2, 3D array
206 @param RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations).
207 @type RT_eigen: numpy rank-2, 3D array
208 @param Ri_prime: The empty rotation matrix for the in-frame rotor motion for state i.
209 @type Ri_prime: numpy rank-2, 3D array
210 @return: The PCS value for the changed position.
211 @rtype: float
212 """
213
214
215 c_sigma = cos(sigma_i)
216 s_sigma = sin(sigma_i)
217 Ri_prime[0, 0] = c_sigma
218 Ri_prime[0, 1] = -s_sigma
219 Ri_prime[1, 0] = s_sigma
220 Ri_prime[1, 1] = c_sigma
221
222
223 R_i = dot(R_eigen, dot(Ri_prime, RT_eigen))
224
225
226 vect = dot(R_i, r_pivot_atom) + r_ln_pivot
227
228
229 length = norm(vect)
230
231
232 proj = dot(vect, dot(A, vect))
233
234
235 pcs = proj / length**5
236
237
238 return pcs
239
240
241 -def pcs_pivot_motion_rotor_qrint(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):
242 """Calculate the PCS value after a pivoted motion for the rotor model.
243
244 @keyword sigma_i: The rotor angle for state i.
245 @type sigma_i: float
246 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor.
247 @type full_in_ref_frame: numpy rank-1 array
248 @keyword r_pivot_atom: The pivot point to atom vector.
249 @type r_pivot_atom: numpy rank-2, 3D array
250 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector.
251 @type r_pivot_atom_rev: numpy rank-2, 3D array
252 @keyword r_ln_pivot: The lanthanide position to pivot point vector.
253 @type r_ln_pivot: numpy rank-2, 3D array
254 @keyword A: The full alignment tensor of the non-moving domain.
255 @type A: numpy rank-2, 3D array
256 @keyword R_eigen: The eigenframe rotation matrix.
257 @type R_eigen: numpy rank-2, 3D array
258 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations).
259 @type RT_eigen: numpy rank-2, 3D array
260 @keyword Ri_prime: The empty rotation matrix for the in-frame rotor motion for state i.
261 @type Ri_prime: numpy rank-2, 3D array
262 @keyword pcs_theta: The storage structure for the back-calculated PCS values.
263 @type pcs_theta: numpy rank-2 array
264 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors.
265 @type pcs_theta_err: numpy rank-2 array
266 @keyword missing_pcs: A structure used to indicate which PCS values are missing.
267 @type missing_pcs: numpy rank-2 array
268 @keyword error_flag: A flag which if True will cause the PCS errors to be estimated and stored in pcs_theta_err.
269 @type error_flag: bool
270 """
271
272
273 c_sigma = cos(sigma_i)
274 s_sigma = sin(sigma_i)
275 Ri_prime[0, 0] = c_sigma
276 Ri_prime[0, 1] = -s_sigma
277 Ri_prime[0, 2] = 0.0
278 Ri_prime[1, 0] = s_sigma
279 Ri_prime[1, 1] = c_sigma
280 Ri_prime[1, 2] = 0.0
281 Ri_prime[2, 0] = 0.0
282 Ri_prime[2, 1] = 0.0
283 Ri_prime[2, 2] = 1.0
284
285
286 R_i = dot(R_eigen, dot(Ri_prime, RT_eigen))
287
288
289 rot_vect_rev = transpose(dot(R_i, r_pivot_atom_rev) + r_ln_pivot)
290 rot_vect = transpose(dot(R_i, r_pivot_atom) + r_ln_pivot)
291
292
293 for j in range(len(r_pivot_atom[0])):
294
295 length_rev = 1.0 / sqrt(inner(rot_vect_rev[j], rot_vect_rev[j]))**5
296 length = 1.0 / sqrt(inner(rot_vect[j], rot_vect[j]))**5
297
298
299 for i in range(len(pcs_theta)):
300
301 if missing_pcs[i, j]:
302 continue
303
304
305 if full_in_ref_frame[i]:
306 proj = dot(rot_vect[j], dot(A[i], rot_vect[j]))
307 length_i = length
308 else:
309 proj = dot(rot_vect_rev[j], dot(A[i], rot_vect_rev[j]))
310 length_i = length_rev
311
312
313 pcs_theta[i, j] += proj * length_i
314
315
316 if error_flag:
317 pcs_theta_err[i, j] += (proj * length_i)**2
318