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