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