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 if dep_check.scipy_module:
31 from scipy.integrate import dblquad
32
33
34 from lib.frame_order.matrix_ops import pcs_pivot_motion_torsionless, pcs_pivot_motion_torsionless_qrint, rotate_daeg
35
36
38 """Generate the rotated 2nd degree Frame Order matrix for the torsionless isotropic cone.
39
40 The cone axis is assumed to be parallel to the z-axis in the eigenframe.
41
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 """
50
51
52 for i in range(9):
53 for j in range(9):
54 matrix[i, j] = 0.0
55
56
57 cos_tmax = cos(cone_theta)
58 cos_tmax2 = cos_tmax**2
59
60
61 matrix[0, 0] = (3.0*cos_tmax2 + 6.0*cos_tmax + 15.0) / 24.0
62 matrix[1, 1] = (cos_tmax2 + 10.0*cos_tmax + 13.0) / 24.0
63 matrix[2, 2] = (4.0*cos_tmax2 + 10.0*cos_tmax + 10.0) / 24.0
64 matrix[3, 3] = matrix[1, 1]
65 matrix[4, 4] = matrix[0, 0]
66 matrix[5, 5] = matrix[2, 2]
67 matrix[6, 6] = matrix[2, 2]
68 matrix[7, 7] = matrix[2, 2]
69 matrix[8, 8] = (cos_tmax2 + cos_tmax + 1.0) / 3.0
70
71
72 matrix[0, 4] = matrix[4, 0] = (cos_tmax2 - 2.0*cos_tmax + 1.0) / 24.0
73 matrix[0, 8] = matrix[8, 0] = -(cos_tmax2 + cos_tmax - 2.0) / 6.0
74 matrix[4, 8] = matrix[8, 4] = matrix[0, 8]
75
76
77 matrix[1, 3] = matrix[3, 1] = matrix[0, 4]
78 matrix[2, 6] = matrix[6, 2] = -matrix[0, 8]
79 matrix[5, 7] = matrix[7, 5] = -matrix[0, 8]
80
81
82 return rotate_daeg(matrix, Rx2_eigen)
83
84
86 """Determine the averaged PCS value via numerical integration.
87
88 @keyword theta_max: The half cone angle.
89 @type theta_max: float
90 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units).
91 @type c: float
92 @keyword r_pivot_atom: The pivot point to atom vector.
93 @type r_pivot_atom: numpy rank-1, 3D array
94 @keyword r_ln_pivot: The lanthanide position to pivot point vector.
95 @type r_ln_pivot: numpy rank-1, 3D array
96 @keyword A: The full alignment tensor of the non-moving domain.
97 @type A: numpy rank-2, 3D array
98 @keyword R_eigen: The eigenframe rotation matrix.
99 @type R_eigen: numpy rank-2, 3D array
100 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations).
101 @type RT_eigen: numpy rank-2, 3D array
102 @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.
103 @type Ri_prime: numpy rank-2, 3D array
104 @return: The averaged PCS value.
105 @rtype: float
106 """
107
108
109 result = dblquad(pcs_pivot_motion_torsionless, -pi, pi, lambda phi: 0.0, lambda phi: theta_max, args=(r_pivot_atom, r_ln_pivot, A, R_eigen, RT_eigen, Ri_prime))
110
111
112 SA = 2.0 * pi * (1.0 - cos(theta_max))
113
114
115 return c * result[0] / SA
116
117
118 -def pcs_numeric_int_iso_cone_torsionless_qrint(points=None, theta_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):
119 """Determine the averaged PCS value via numerical integration.
120
121 @keyword points: The Sobol points in the torsion-tilt angle space.
122 @type points: numpy rank-2, 3D array
123 @keyword theta_max: The half cone angle.
124 @type theta_max: float
125 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units).
126 @type c: numpy rank-1 array
127 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor.
128 @type full_in_ref_frame: numpy rank-1 array
129 @keyword r_pivot_atom: The pivot point to atom vector.
130 @type r_pivot_atom: numpy rank-2, 3D array
131 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector.
132 @type r_pivot_atom_rev: numpy rank-2, 3D array
133 @keyword r_ln_pivot: The lanthanide position to pivot point vector.
134 @type r_ln_pivot: numpy rank-2, 3D array
135 @keyword A: The full alignment tensor of the non-moving domain.
136 @type A: numpy rank-2, 3D array
137 @keyword R_eigen: The eigenframe rotation matrix.
138 @type R_eigen: numpy rank-2, 3D array
139 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations).
140 @type RT_eigen: numpy rank-2, 3D array
141 @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.
142 @type Ri_prime: numpy rank-2, 3D array
143 @keyword pcs_theta: The storage structure for the back-calculated PCS values.
144 @type pcs_theta: numpy rank-2 array
145 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors.
146 @type pcs_theta_err: numpy rank-2 array
147 @keyword missing_pcs: A structure used to indicate which PCS values are missing.
148 @type missing_pcs: numpy rank-2 array
149 @keyword error_flag: A flag which if True will cause the PCS errors to be estimated and stored in pcs_theta_err.
150 @type error_flag: bool
151 """
152
153
154 for i in range(len(pcs_theta)):
155 for j in range(len(pcs_theta[i])):
156 pcs_theta[i, j] = 0.0
157 pcs_theta_err[i, j] = 0.0
158
159
160 num = 0
161 for i in range(len(points)):
162
163 theta, phi = points[i]
164
165
166 if theta > theta_max:
167 continue
168
169
170 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)
171
172
173 num += 1
174
175
176 for i in range(len(pcs_theta)):
177 for j in range(len(pcs_theta[i])):
178
179 pcs_theta[i, j] = c[i] * pcs_theta[i, j] / float(num)
180
181
182 if error_flag:
183 pcs_theta_err[i, j] = abs(pcs_theta_err[i, j] / float(num) - pcs_theta[i, j]**2) / float(num)
184 pcs_theta_err[i, j] = c[i] * sqrt(pcs_theta_err[i, j])
185 print("%8.3f +/- %-8.3f" % (pcs_theta[i, j]*1e6, pcs_theta_err[i, j]*1e6))
186