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, sqrt
30 from numpy import sinc
31 if dep_check.scipy_module:
32 from scipy.integrate import tplquad
33
34
35 from lib.frame_order.matrix_ops import pcs_pivot_motion_full, pcs_pivot_motion_full_qrint, rotate_daeg
36
37
39 """Generate the rotated 2nd degree Frame Order matrix for the isotropic cone.
40
41 The cone axis is assumed to be parallel to the z-axis in the eigenframe.
42
43 @param matrix: The Frame Order matrix, 2nd degree to be populated.
44 @type matrix: numpy 9D, rank-2 array
45 @param Rx2_eigen: The Kronecker product of the eigenframe rotation matrix with itself.
46 @type Rx2_eigen: numpy 9D, rank-2 array
47 @param cone_theta: The cone opening angle.
48 @type cone_theta: float
49 @param sigma_max: The maximum torsion angle.
50 @type sigma_max: float
51 """
52
53
54 populate_2nd_eigenframe_iso_cone(matrix, cone_theta, sigma_max)
55
56
57 return rotate_daeg(matrix, Rx2_eigen)
58
59
60 -def pcs_numeric_int_iso_cone(theta_max=None, sigma_max=None, c=None, r_pivot_atom=None, r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None):
61 """Determine the averaged PCS value via numerical integration.
62
63 @keyword theta_max: The half cone angle.
64 @type theta_max: float
65 @keyword sigma_max: The maximum torsion angle.
66 @type sigma_max: float
67 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units).
68 @type c: float
69 @keyword r_pivot_atom: The pivot point to atom vector.
70 @type r_pivot_atom: numpy rank-1, 3D array
71 @keyword r_ln_pivot: The lanthanide position to pivot point vector.
72 @type r_ln_pivot: numpy rank-1, 3D array
73 @keyword A: The full alignment tensor of the non-moving domain.
74 @type A: numpy rank-2, 3D array
75 @keyword R_eigen: The eigenframe rotation matrix.
76 @type R_eigen: numpy rank-2, 3D array
77 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations).
78 @type RT_eigen: numpy rank-2, 3D array
79 @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.
80 @type Ri_prime: numpy rank-2, 3D array
81 @return: The averaged PCS value.
82 @rtype: float
83 """
84
85
86 result = tplquad(pcs_pivot_motion_full, -sigma_max, sigma_max, lambda phi: -pi, lambda phi: pi, lambda theta, phi: 0.0, lambda theta, phi: theta_max, args=(r_pivot_atom, r_ln_pivot, A, R_eigen, RT_eigen, Ri_prime))
87
88
89 SA = 4.0 * pi * sigma_max * (1.0 - cos(theta_max))
90
91
92 return c * result[0] / SA
93
94
95 -def pcs_numeric_int_iso_cone_qrint(points=None, theta_max=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):
96 """Determine the averaged PCS value via numerical integration.
97
98 @keyword points: The Sobol points in the torsion-tilt angle space.
99 @type points: numpy rank-2, 3D array
100 @keyword theta_max: The half cone angle.
101 @type theta_max: float
102 @keyword sigma_max: The maximum torsion angle.
103 @type sigma_max: float
104 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units).
105 @type c: numpy rank-1 array
106 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor.
107 @type full_in_ref_frame: numpy rank-1 array
108 @keyword r_pivot_atom: The pivot point to atom vector.
109 @type r_pivot_atom: numpy rank-2, 3D array
110 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector.
111 @type r_pivot_atom_rev: numpy rank-2, 3D array
112 @keyword r_ln_pivot: The lanthanide position to pivot point vector.
113 @type r_ln_pivot: numpy rank-2, 3D array
114 @keyword A: The full alignment tensor of the non-moving domain.
115 @type A: numpy rank-2, 3D array
116 @keyword R_eigen: The eigenframe rotation matrix.
117 @type R_eigen: numpy rank-2, 3D array
118 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations).
119 @type RT_eigen: numpy rank-2, 3D array
120 @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.
121 @type Ri_prime: numpy rank-2, 3D array
122 @keyword pcs_theta: The storage structure for the back-calculated PCS values.
123 @type pcs_theta: numpy rank-2 array
124 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors.
125 @type pcs_theta_err: numpy rank-2 array
126 @keyword missing_pcs: A structure used to indicate which PCS values are missing.
127 @type missing_pcs: numpy rank-2 array
128 @keyword error_flag: A flag which if True will cause the PCS errors to be estimated and stored in pcs_theta_err.
129 @type error_flag: bool
130 """
131
132
133 for i in range(len(pcs_theta)):
134 for j in range(len(pcs_theta[i])):
135 pcs_theta[i, j] = 0.0
136 pcs_theta_err[i, j] = 0.0
137
138
139 num = 0
140 for i in range(len(points)):
141
142 theta, phi, sigma = points[i]
143
144
145 if theta > theta_max:
146 continue
147 if sigma > sigma_max or sigma < -sigma_max:
148 continue
149
150
151 pcs_pivot_motion_full_qrint(theta_i=theta, phi_i=phi, 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)
152
153
154 num += 1
155
156
157 for i in range(len(pcs_theta)):
158 for j in range(len(pcs_theta[i])):
159
160 pcs_theta[i, j] = c[i] * pcs_theta[i, j] / float(num)
161
162
163 if error_flag:
164 pcs_theta_err[i, j] = abs(pcs_theta_err[i, j] / float(num) - pcs_theta[i, j]**2) / float(num)
165 pcs_theta_err[i, j] = c[i] * sqrt(pcs_theta_err[i, j])
166 print("%8.3f +/- %-8.3f" % (pcs_theta[i, j]*1e6, pcs_theta_err[i, j]*1e6))
167
168
170 """Populate the 1st degree Frame Order matrix in the eigenframe for an isotropic cone.
171
172 The cone axis is assumed to be parallel to the z-axis in the eigenframe.
173
174 @param matrix: The Frame Order matrix, 1st degree.
175 @type matrix: numpy 3D, rank-2 array
176 @param angle: The cone angle.
177 @type angle: float
178 """
179
180
181 for i in range(3):
182 for j in range(3):
183 matrix[i, j] = 0.0
184
185
186 matrix[2, 2] = (cos(angle) + 1.0) / 2.0
187
188
190 """Populate the 2nd degree Frame Order matrix in the eigenframe for the isotropic cone.
191
192 The cone axis is assumed to be parallel to the z-axis in the eigenframe.
193
194
195 @param matrix: The Frame Order matrix, 2nd degree.
196 @type matrix: numpy 9D, rank-2 array
197 @param tmax: The cone opening angle.
198 @type tmax: float
199 @param smax: The maximum torsion angle.
200 @type smax: float
201 """
202
203
204 for i in range(9):
205 for j in range(9):
206 matrix[i, j] = 0.0
207
208
209 sinc_smax = sinc(smax/pi)
210 sinc_2smax = sinc(2.0*smax/pi)
211 cos_tmax = cos(tmax)
212 cos_tmax2 = cos_tmax**2
213
214
215 fact_sinc_2smax = sinc_2smax*(cos_tmax2 + 4.0*cos_tmax + 7.0) / 24.0
216 fact_cos_tmax2 = (cos_tmax2 + cos_tmax + 4.0) / 12.0
217 fact_cos_tmax = (cos_tmax + 1.0) / 4.0
218
219
220 matrix[0, 0] = fact_sinc_2smax + fact_cos_tmax2
221 matrix[1, 1] = fact_sinc_2smax + fact_cos_tmax
222 matrix[2, 2] = sinc_smax * (2.0*cos_tmax2 + 5.0*cos_tmax + 5.0) / 12.0
223 matrix[3, 3] = matrix[1, 1]
224 matrix[4, 4] = matrix[0, 0]
225 matrix[5, 5] = matrix[2, 2]
226 matrix[6, 6] = matrix[2, 2]
227 matrix[7, 7] = matrix[2, 2]
228 matrix[8, 8] = (cos_tmax2 + cos_tmax + 1.0) / 3.0
229
230
231 matrix[0, 4] = matrix[4, 0] = -fact_sinc_2smax + fact_cos_tmax2
232 matrix[0, 8] = matrix[8, 0] = -(cos_tmax2 + cos_tmax - 2.0) / 6.0
233 matrix[4, 8] = matrix[8, 4] = matrix[0, 8]
234
235
236 matrix[1, 3] = matrix[3, 1] = fact_sinc_2smax - fact_cos_tmax
237 matrix[2, 6] = matrix[6, 2] = sinc_smax * (cos_tmax2 + cos_tmax - 2.0) / 6.0
238 matrix[5, 7] = matrix[7, 5] = matrix[2, 6]
239