lib.frame_order.pseudo_ellipse

1 ############################################################################### 2 # # 3 # Copyright (C) 2009-2013 Edward d'Auvergne # 4 # # 5 # This file is part of the program relax (http://www.nmr-relax.com). # 6 # # 7 # This program is free software: you can redistribute it and/or modify # 8 # it under the terms of the GNU General Public License as published by # 9 # the Free Software Foundation, either version 3 of the License, or # 10 # (at your option) any later version. # 11 # # 12 # This program is distributed in the hope that it will be useful, # 13 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 15 # GNU General Public License for more details. # 16 # # 17 # You should have received a copy of the GNU General Public License # 18 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 19 # # 20 ############################################################################### 21 22 # Module docstring. 23 """Module for the handling of Frame Order.""" 24 25 # Dependency check module. 26 import dep_check 27 28 # Python module imports. 29 from math import cos, pi, sin, sqrt 30 if dep_check.scipy_module: 31 from scipy.integrate import dblquad, quad 32 33 # relax module imports. 34 from lib.geometry.pec import pec 35 from lib.frame_order.matrix_ops import pcs_pivot_motion_torsionless, pcs_pivot_motion_torsionless_qrint, rotate_daeg 36 from lib.frame_order.pseudo_ellipse import tmax_pseudo_ellipse 37 38

39 -def compile_2nd_matrix_pseudo_ellipse_torsionless(matrix, Rx2_eigen, theta_x, theta_y):

40 """Generate the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 41 42 @param matrix: The Frame Order matrix, 2nd degree to be populated. 43 @type matrix: numpy 9D, rank-2 array 44 @param Rx2_eigen: The Kronecker product of the eigenframe rotation matrix with itself. 45 @type Rx2_eigen: numpy 9D, rank-2 array 46 @param theta_x: The cone opening angle along x. 47 @type theta_x: float 48 @param theta_y: The cone opening angle along y. 49 @type theta_y: float 50 """ 51 52 # The surface area normalisation factor. 53 fact = 1.0 / (6.0 * pec(theta_x, theta_y)) 54 55 # Diagonal. 56 matrix[0, 0] = fact * (6.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_00, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 57 matrix[1, 1] = fact * (2.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_11, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 58 matrix[2, 2] = fact * (5.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_22, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 59 matrix[3, 3] = matrix[1, 1] 60 matrix[4, 4] = fact * (6.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_44, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 61 matrix[5, 5] = fact * (5.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_55, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 62 matrix[6, 6] = matrix[2, 2] 63 matrix[7, 7] = matrix[5, 5] 64 matrix[8, 8] = fact * quad(part_int_daeg2_pseudo_ellipse_torsionless_88, -pi, pi, args=(theta_x, theta_y), full_output=1)[0] 65 66 # Off diagonal set 1. 67 matrix[0, 4] = matrix[4, 0] = fact * (2.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_04, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 68 matrix[0, 8] = matrix[8, 0] = fact * (4.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_08, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 69 matrix[4, 8] = matrix[8, 4] = fact * (4.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_48, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 70 71 # Off diagonal set 2. 72 matrix[1, 3] = matrix[3, 1] = matrix[0, 4] 73 matrix[2, 6] = matrix[6, 2] = -matrix[0, 8] 74 matrix[5, 7] = matrix[7, 5] = -matrix[4, 8] 75 76 # Rotate and return the frame order matrix. 77 return rotate_daeg(matrix, Rx2_eigen)

78 79

80 -def part_int_daeg2_pseudo_ellipse_torsionless_00(phi, x, y):

81 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 82 83 @param phi: The azimuthal tilt-torsion angle. 84 @type phi: float 85 @param x: The cone opening angle along x. 86 @type x: float 87 @param y: The cone opening angle along y. 88 @type y: float 89 @return: The theta partial integral. 90 @rtype: float 91 """ 92 93 # Theta max. 94 tmax = tmax_pseudo_ellipse(phi, x, y) 95 96 # The theta integral. 97 return (2*cos(phi)**4*cos(tmax) + 6*cos(phi)**2*sin(phi)**2)*sin(tmax)**2 - (6*sin(phi)**4 + 2*cos(phi)**4)*cos(tmax)

98 99

100 -def part_int_daeg2_pseudo_ellipse_torsionless_04(phi, x, y):

101 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 102 103 @param phi: The azimuthal tilt-torsion angle. 104 @type phi: float 105 @param x: The cone opening angle along x. 106 @type x: float 107 @param y: The cone opening angle along y. 108 @type y: float 109 @return: The theta partial integral. 110 @rtype: float 111 """ 112 113 # Theta max. 114 tmax = tmax_pseudo_ellipse(phi, x, y) 115 116 # The theta integral. 117 return (2*cos(phi)**2*sin(phi)**2*cos(tmax) - 6*cos(phi)**2*sin(phi)**2)*sin(tmax)**2 - 8*cos(phi)**2*sin(phi)**2*cos(tmax)

118 119

120 -def part_int_daeg2_pseudo_ellipse_torsionless_08(phi, x, y):

121 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 122 123 @param phi: The azimuthal tilt-torsion angle. 124 @type phi: float 125 @param x: The cone opening angle along x. 126 @type x: float 127 @param y: The cone opening angle along y. 128 @type y: float 129 @return: The theta partial integral. 130 @rtype: float 131 """ 132 133 # Theta max. 134 tmax = tmax_pseudo_ellipse(phi, x, y) 135 136 # The theta integral. 137 return 2*cos(phi)**2*cos(tmax)**3 - 6*cos(phi)**2*cos(tmax)

138 139

140 -def part_int_daeg2_pseudo_ellipse_torsionless_11(phi, x, y):

141 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 142 143 @param phi: The azimuthal tilt-torsion angle. 144 @type phi: float 145 @param x: The cone opening angle along x. 146 @type x: float 147 @param y: The cone opening angle along y. 148 @type y: float 149 @return: The theta partial integral. 150 @rtype: float 151 """ 152 153 # Theta max. 154 tmax = tmax_pseudo_ellipse(phi, x, y) 155 156 # The theta integral. 157 return (2*cos(phi)**2*sin(phi)**2*cos(tmax) + 3*sin(phi)**4 + 3*cos(phi)**4)*sin(tmax)**2 - 8*cos(phi)**2*sin(phi)**2*cos(tmax)

158 159

160 -def part_int_daeg2_pseudo_ellipse_torsionless_22(phi, x, y):

161 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 162 163 @param phi: The azimuthal tilt-torsion angle. 164 @type phi: float 165 @param x: The cone opening angle along x. 166 @type x: float 167 @param y: The cone opening angle along y. 168 @type y: float 169 @return: The theta partial integral. 170 @rtype: float 171 """ 172 173 # Theta max. 174 tmax = tmax_pseudo_ellipse(phi, x, y) 175 176 # The theta integral. 177 return (2*sin(phi)**2 - 2)*cos(tmax)**3 - 3*sin(phi)**2*cos(tmax)**2

178 179

180 -def part_int_daeg2_pseudo_ellipse_torsionless_44(phi, x, y):

181 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 182 183 @param phi: The azimuthal tilt-torsion angle. 184 @type phi: float 185 @param x: The cone opening angle along x. 186 @type x: float 187 @param y: The cone opening angle along y. 188 @type y: float 189 @return: The theta partial integral. 190 @rtype: float 191 """ 192 193 # Theta max. 194 tmax = tmax_pseudo_ellipse(phi, x, y) 195 196 # The theta integral. 197 return (2*sin(phi)**4*cos(tmax) + 6*cos(phi)**2*sin(phi)**2)*sin(tmax)**2 - (2*sin(phi)**4 + 6*cos(phi)**4)*cos(tmax)

198 199

200 -def part_int_daeg2_pseudo_ellipse_torsionless_48(phi, x, y):

201 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 202 203 @param phi: The azimuthal tilt-torsion angle. 204 @type phi: float 205 @param x: The cone opening angle along x. 206 @type x: float 207 @param y: The cone opening angle along y. 208 @type y: float 209 @return: The theta partial integral. 210 @rtype: float 211 """ 212 213 # Theta max. 214 tmax = tmax_pseudo_ellipse(phi, x, y) 215 216 # The theta integral. 217 return 2*sin(phi)**2*cos(tmax)**3 - 6*sin(phi)**2*cos(tmax)

218 219

220 -def part_int_daeg2_pseudo_ellipse_torsionless_55(phi, x, y):

221 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 222 223 @param phi: The azimuthal tilt-torsion angle. 224 @type phi: float 225 @param x: The cone opening angle along x. 226 @type x: float 227 @param y: The cone opening angle along y. 228 @type y: float 229 @return: The theta partial integral. 230 @rtype: float 231 """ 232 233 # Theta max. 234 tmax = tmax_pseudo_ellipse(phi, x, y) 235 236 # The theta integral. 237 return (2*cos(phi)**2 - 2)*cos(tmax)**3 - 3*cos(phi)**2*cos(tmax)**2

238 239

240 -def part_int_daeg2_pseudo_ellipse_torsionless_88(phi, x, y):

241 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 242 243 @param phi: The azimuthal tilt-torsion angle. 244 @type phi: float 245 @param x: The cone opening angle along x. 246 @type x: float 247 @param y: The cone opening angle along y. 248 @type y: float 249 @return: The theta partial integral. 250 @rtype: float 251 """ 252 253 # Theta max. 254 tmax = tmax_pseudo_ellipse(phi, x, y) 255 256 # The theta integral. 257 return 2 - 2*cos(tmax)**3

258 259

260 -def pcs_numeric_int_pseudo_ellipse_torsionless(theta_x=None, theta_y=None, c=None, r_pivot_atom=None, r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None):

261 """Determine the averaged PCS value via numerical integration. 262 263 @keyword theta_x: The x-axis half cone angle. 264 @type theta_x: float 265 @keyword theta_y: The y-axis half cone angle. 266 @type theta_y: float 267 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units). 268 @type c: float 269 @keyword r_pivot_atom: The pivot point to atom vector. 270 @type r_pivot_atom: numpy rank-1, 3D array 271 @keyword r_ln_pivot: The lanthanide position to pivot point vector. 272 @type r_ln_pivot: numpy rank-1, 3D array 273 @keyword A: The full alignment tensor of the non-moving domain. 274 @type A: numpy rank-2, 3D array 275 @keyword R_eigen: The eigenframe rotation matrix. 276 @type R_eigen: numpy rank-2, 3D array 277 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations). 278 @type RT_eigen: numpy rank-2, 3D array 279 @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. 280 @type Ri_prime: numpy rank-2, 3D array 281 @return: The averaged PCS value. 282 @rtype: float 283 """ 284 285 def pseudo_ellipse(phi): 286 """The pseudo-ellipse wrapper formula.""" 287 288 return tmax_pseudo_ellipse(phi, theta_x, theta_y)

289 290 # Perform numerical integration. 291 result = dblquad(pcs_pivot_motion_torsionless, -pi, pi, lambda phi: 0.0, pseudo_ellipse, args=(r_pivot_atom, r_ln_pivot, A, R_eigen, RT_eigen, Ri_prime)) 292 293 # The surface area normalisation factor. 294 SA = pec(theta_x, theta_y) 295 296 # Return the value. 297 return c * result[0] / SA 298 299

300 -def pcs_numeric_int_pseudo_ellipse_torsionless_qrint(points=None, theta_x=None, theta_y=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):

301 """Determine the averaged PCS value via numerical integration. 302 303 @keyword points: The Sobol points in the torsion-tilt angle space. 304 @type points: numpy rank-2, 3D array 305 @keyword theta_x: The x-axis half cone angle. 306 @type theta_x: float 307 @keyword theta_y: The y-axis half cone angle. 308 @type theta_y: float 309 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units). 310 @type c: numpy rank-1 array 311 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 312 @type full_in_ref_frame: numpy rank-1 array 313 @keyword r_pivot_atom: The pivot point to atom vector. 314 @type r_pivot_atom: numpy rank-2, 3D array 315 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector. 316 @type r_pivot_atom_rev: numpy rank-2, 3D array 317 @keyword r_ln_pivot: The lanthanide position to pivot point vector. 318 @type r_ln_pivot: numpy rank-2, 3D array 319 @keyword A: The full alignment tensor of the non-moving domain. 320 @type A: numpy rank-2, 3D array 321 @keyword R_eigen: The eigenframe rotation matrix. 322 @type R_eigen: numpy rank-2, 3D array 323 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations). 324 @type RT_eigen: numpy rank-2, 3D array 325 @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. 326 @type Ri_prime: numpy rank-2, 3D array 327 @keyword pcs_theta: The storage structure for the back-calculated PCS values. 328 @type pcs_theta: numpy rank-2 array 329 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors. 330 @type pcs_theta_err: numpy rank-2 array 331 @keyword missing_pcs: A structure used to indicate which PCS values are missing. 332 @type missing_pcs: numpy rank-2 array 333 @keyword error_flag: A flag which if True will cause the PCS errors to be estimated and stored in pcs_theta_err. 334 @type error_flag: bool 335 """ 336 337 # Clear the data structures. 338 for i in range(len(pcs_theta)): 339 for j in range(len(pcs_theta[i])): 340 pcs_theta[i, j] = 0.0 341 pcs_theta_err[i, j] = 0.0 342 343 # Loop over the samples. 344 num = 0 345 for i in range(len(points)): 346 # Unpack the point. 347 theta, phi = points[i] 348 349 # Calculate theta_max. 350 theta_max = tmax_pseudo_ellipse(phi, theta_x, theta_y) 351 352 # Outside of the distribution, so skip the point. 353 if theta > theta_max: 354 continue 355 356 # Calculate the PCSs for this state. 357 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) 358 359 # Increment the number of points. 360 num += 1 361 362 # Calculate the PCS and error. 363 for i in range(len(pcs_theta)): 364 for j in range(len(pcs_theta[i])): 365 # The average PCS. 366 pcs_theta[i, j] = c[i] * pcs_theta[i, j] / float(num) 367 368 # The error. 369 if error_flag: 370 pcs_theta_err[i, j] = abs(pcs_theta_err[i, j] / float(num) - pcs_theta[i, j]**2) / float(num) 371 pcs_theta_err[i, j] = c[i] * sqrt(pcs_theta_err[i, j]) 372 print("%8.3f +/- %-8.3f" % (pcs_theta[i, j]*1e6, pcs_theta_err[i, j]*1e6))

373

Source Code for Module lib.frame_order.pseudo_ellipse_torsionless