lib.frame_order.pseudo_ellipse

1 ############################################################################### 2 # # 3 # Copyright (C) 2009-2012,2014 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 # Python module imports. 26 from math import cos, pi, sin 27 from numpy import divide, dot, eye, float64, multiply, swapaxes, tensordot 28 try: 29 from scipy.integrate import dblquad, quad 30 except ImportError: 31 pass 32 33 # relax module imports. 34 from lib.geometry.pec import pec 35 from lib.frame_order.matrix_ops import pcs_pivot_motion_torsionless_qr_int, pcs_pivot_motion_torsionless_quad_int, rotate_daeg 36 from lib.frame_order.pseudo_ellipse import tmax_pseudo_ellipse, tmax_pseudo_ellipse_array 37 38

39 -def compile_1st_matrix_pseudo_ellipse_torsionless(matrix, R_eigen, theta_x, theta_y):

40 """Generate the 1st degree Frame Order matrix for the torsionless pseudo-ellipse. 41 42 @param matrix: The Frame Order matrix, 1st degree to be populated. 43 @type matrix: numpy 3D, rank-2 array 44 @param R_eigen: The eigenframe rotation matrix. 45 @type R_eigen: numpy 3D, 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 = 2.0 * pec(theta_x, theta_y) 54 55 # Invert. 56 if fact == 0.0: 57 fact = 1e100 58 else: 59 fact = 1.0 / fact 60 61 # Numerical integration of phi of each element. 62 matrix[0, 0] = fact * (2.0*pi + quad(part_int_daeg1_pseudo_ellipse_00, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 63 matrix[1, 1] = fact * (2.0*pi + quad(part_int_daeg1_pseudo_ellipse_11, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 64 matrix[2, 2] = fact * quad(part_int_daeg1_pseudo_ellipse_22, -pi, pi, args=(theta_x, theta_y), full_output=1)[0] 65 66 # Rotate and return the frame order matrix. 67 return rotate_daeg(matrix, R_eigen)

68 69

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

71 """Generate the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 72 73 @param matrix: The Frame Order matrix, 2nd degree to be populated. 74 @type matrix: numpy 9D, rank-2 array 75 @param Rx2_eigen: The Kronecker product of the eigenframe rotation matrix with itself. 76 @type Rx2_eigen: numpy 9D, rank-2 array 77 @param theta_x: The cone opening angle along x. 78 @type theta_x: float 79 @param theta_y: The cone opening angle along y. 80 @type theta_y: float 81 """ 82 83 # The rigid case. 84 if theta_x == 0.0: 85 # Set up the matrix as the identity. 86 matrix[:] = 0.0 87 for i in range(len(matrix)): 88 matrix[i, i] = 1.0 89 90 # Rotate and return the frame order matrix. 91 return rotate_daeg(matrix, Rx2_eigen) 92 93 # The surface area normalisation factor. 94 fact = 6.0 * pec(theta_x, theta_y) 95 fact2 = 0.5 * fact 96 97 # Invert. 98 if fact == 0.0: 99 fact = 1e100 100 fact2 = 1e100 101 else: 102 fact = 1.0 / fact 103 fact2 = 1.0 / fact2 104 105 # Diagonal. 106 matrix[0, 0] = fact2 * (3.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_00, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 107 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]) 108 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]) 109 matrix[3, 3] = matrix[1, 1] 110 matrix[4, 4] = fact2 * (3.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_44, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 111 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]) 112 matrix[6, 6] = matrix[2, 2] 113 matrix[7, 7] = matrix[5, 5] 114 matrix[8, 8] = fact2 * (2.0*pi - quad(part_int_daeg2_pseudo_ellipse_torsionless_88, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 115 116 # Off diagonal set 1. 117 matrix[0, 4] = matrix[4, 0] = fact2 * (pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_04, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 118 matrix[0, 8] = matrix[8, 0] = fact2 * (2.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_08, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 119 matrix[4, 8] = matrix[8, 4] = fact2 * (2.0*pi + quad(part_int_daeg2_pseudo_ellipse_torsionless_48, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 120 121 # Off diagonal set 2. 122 matrix[1, 3] = matrix[3, 1] = matrix[0, 4] 123 matrix[2, 6] = matrix[6, 2] = -matrix[0, 8] 124 matrix[5, 7] = matrix[7, 5] = -matrix[4, 8] 125 126 # Rotate and return the frame order matrix. 127 return rotate_daeg(matrix, Rx2_eigen)

128 129

130 -def part_int_daeg1_pseudo_ellipse_00(phi, x, y):

131 """The theta-sigma partial integral of the 1st degree Frame Order matrix element 00 for the pseudo-ellipse. 132 133 @param phi: The azimuthal tilt-torsion angle. 134 @type phi: float 135 @param x: The cone opening angle along x. 136 @type x: float 137 @param y: The cone opening angle along y. 138 @type y: float 139 @return: The theta-sigma partial integral. 140 @rtype: float 141 """ 142 143 # Theta max. 144 tmax = tmax_pseudo_ellipse(phi, x, y) 145 146 # The theta-sigma integral. 147 return cos(phi)**2 * sin(tmax)**2 - 2.0 * sin(phi)**2 * cos(tmax)

148 149

150 -def part_int_daeg1_pseudo_ellipse_11(phi, x, y):

151 """The theta-sigma partial integral of the 1st degree Frame Order matrix element 11 for the pseudo-ellipse. 152 153 @param phi: The azimuthal tilt-torsion angle. 154 @type phi: float 155 @param x: The cone opening angle along x. 156 @type x: float 157 @param y: The cone opening angle along y. 158 @type y: float 159 @return: The theta-sigma partial integral. 160 @rtype: float 161 """ 162 163 # Theta max. 164 tmax = tmax_pseudo_ellipse(phi, x, y) 165 166 # The theta-sigma integral. 167 return sin(phi)**2 * sin(tmax)**2 - 2.0 * cos(phi)**2 * cos(tmax)

168 169

170 -def part_int_daeg1_pseudo_ellipse_22(phi, x, y):

171 """The theta-sigma partial integral of the 1st degree Frame Order matrix element 22 for the pseudo-ellipse. 172 173 @param phi: The azimuthal tilt-torsion angle. 174 @type phi: float 175 @param x: The cone opening angle along x. 176 @type x: float 177 @param y: The cone opening angle along y. 178 @type y: float 179 @return: The theta-sigma partial integral. 180 @rtype: float 181 """ 182 183 # Theta max. 184 tmax = tmax_pseudo_ellipse(phi, x, y) 185 186 # The theta-sigma integral. 187 return sin(tmax)**2

188 189

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

191 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 192 193 @param phi: The azimuthal tilt-torsion angle. 194 @type phi: float 195 @param x: The cone opening angle along x. 196 @type x: float 197 @param y: The cone opening angle along y. 198 @type y: float 199 @return: The theta partial integral. 200 @rtype: float 201 """ 202 203 # Theta max. 204 tmax = tmax_pseudo_ellipse(phi, x, y) 205 206 # The theta integral. 207 return (cos(phi)**4*cos(tmax) + 3.0*cos(phi)**2.0*sin(phi)**2)*sin(tmax)**2 - (3.0*sin(phi)**4 + cos(phi)**4)*cos(tmax)

208 209

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

211 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 212 213 @param phi: The azimuthal tilt-torsion angle. 214 @type phi: float 215 @param x: The cone opening angle along x. 216 @type x: float 217 @param y: The cone opening angle along y. 218 @type y: float 219 @return: The theta partial integral. 220 @rtype: float 221 """ 222 223 # Theta max. 224 tmax = tmax_pseudo_ellipse(phi, x, y) 225 226 # The theta integral. 227 return (cos(phi)**2*sin(phi)**2*cos(tmax) - 3.0*cos(phi)**2*sin(phi)**2)*sin(tmax)**2 - 4.0*cos(phi)**2*sin(phi)**2*cos(tmax)

228 229

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

231 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 232 233 @param phi: The azimuthal tilt-torsion angle. 234 @type phi: float 235 @param x: The cone opening angle along x. 236 @type x: float 237 @param y: The cone opening angle along y. 238 @type y: float 239 @return: The theta partial integral. 240 @rtype: float 241 """ 242 243 # Theta max. 244 tmax = tmax_pseudo_ellipse(phi, x, y) 245 246 # The theta integral. 247 return cos(phi)**2*cos(tmax)**3 - 3.0*cos(phi)**2*cos(tmax)

248 249

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

251 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 252 253 @param phi: The azimuthal tilt-torsion angle. 254 @type phi: float 255 @param x: The cone opening angle along x. 256 @type x: float 257 @param y: The cone opening angle along y. 258 @type y: float 259 @return: The theta partial integral. 260 @rtype: float 261 """ 262 263 # Theta max. 264 tmax = tmax_pseudo_ellipse(phi, x, y) 265 266 # The theta integral. 267 return (2.0*cos(phi)**2*sin(phi)**2*cos(tmax) + 3.0*sin(phi)**4 + 3.0*cos(phi)**4)*sin(tmax)**2 - 8.0*cos(phi)**2*sin(phi)**2*cos(tmax)

268 269

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

271 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 272 273 @param phi: The azimuthal tilt-torsion angle. 274 @type phi: float 275 @param x: The cone opening angle along x. 276 @type x: float 277 @param y: The cone opening angle along y. 278 @type y: float 279 @return: The theta partial integral. 280 @rtype: float 281 """ 282 283 # Theta max. 284 tmax = tmax_pseudo_ellipse(phi, x, y) 285 286 # The theta integral. 287 return 2.0*cos(phi)**2*cos(tmax)**3 + 3.0*sin(phi)**2*cos(tmax)**2

288 289

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

291 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 292 293 @param phi: The azimuthal tilt-torsion angle. 294 @type phi: float 295 @param x: The cone opening angle along x. 296 @type x: float 297 @param y: The cone opening angle along y. 298 @type y: float 299 @return: The theta partial integral. 300 @rtype: float 301 """ 302 303 # Theta max. 304 tmax = tmax_pseudo_ellipse(phi, x, y) 305 306 # The theta integral. 307 return (sin(phi)**4*cos(tmax) + 3.0*cos(phi)**2*sin(phi)**2)*sin(tmax)**2 - (sin(phi)**4 + 3.0*cos(phi)**4)*cos(tmax)

308 309

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

311 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 312 313 @param phi: The azimuthal tilt-torsion angle. 314 @type phi: float 315 @param x: The cone opening angle along x. 316 @type x: float 317 @param y: The cone opening angle along y. 318 @type y: float 319 @return: The theta partial integral. 320 @rtype: float 321 """ 322 323 # Theta max. 324 tmax = tmax_pseudo_ellipse(phi, x, y) 325 326 # The theta integral. 327 return sin(phi)**2*cos(tmax)**3 - 3.0*sin(phi)**2*cos(tmax)

328 329

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

331 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 332 333 @param phi: The azimuthal tilt-torsion angle. 334 @type phi: float 335 @param x: The cone opening angle along x. 336 @type x: float 337 @param y: The cone opening angle along y. 338 @type y: float 339 @return: The theta partial integral. 340 @rtype: float 341 """ 342 343 # Theta max. 344 tmax = tmax_pseudo_ellipse(phi, x, y) 345 346 # The theta integral. 347 return 2.0*sin(phi)**2*cos(tmax)**3 + 3.0*cos(phi)**2*cos(tmax)**2

348 349

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

351 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 352 353 @param phi: The azimuthal tilt-torsion angle. 354 @type phi: float 355 @param x: The cone opening angle along x. 356 @type x: float 357 @param y: The cone opening angle along y. 358 @type y: float 359 @return: The theta partial integral. 360 @rtype: float 361 """ 362 363 # Theta max. 364 tmax = tmax_pseudo_ellipse(phi, x, y) 365 366 # The theta integral. 367 return cos(tmax)**3

368 369

370 -def pcs_numeric_qr_int_pseudo_ellipse_torsionless(points=None, max_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):

371 """Determine the averaged PCS value via numerical integration. 372 373 @keyword points: The Sobol points in the torsion-tilt angle space. 374 @type points: numpy rank-2, 3D array 375 @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. 376 @type max_points: int 377 @keyword theta_x: The x-axis half cone angle. 378 @type theta_x: float 379 @keyword theta_y: The y-axis half cone angle. 380 @type theta_y: float 381 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units). 382 @type c: numpy rank-1 array 383 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 384 @type full_in_ref_frame: numpy rank-1 array 385 @keyword r_pivot_atom: The pivot point to atom vector. 386 @type r_pivot_atom: numpy rank-2, 3D array 387 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector. 388 @type r_pivot_atom_rev: numpy rank-2, 3D array 389 @keyword r_ln_pivot: The lanthanide position to pivot point vector. 390 @type r_ln_pivot: numpy rank-2, 3D array 391 @keyword A: The full alignment tensor of the non-moving domain. 392 @type A: numpy rank-2, 3D array 393 @keyword R_eigen: The eigenframe rotation matrix. 394 @type R_eigen: numpy rank-2, 3D array 395 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations). 396 @type RT_eigen: numpy rank-2, 3D array 397 @keyword Ri_prime: The array of pre-calculated rotation matrices for the in-frame torsionless, pseudo-elliptic cone motion, used to calculate the PCS for each state i in the numerical integration. 398 @type Ri_prime: numpy rank-3, array of 3D arrays 399 @keyword pcs_theta: The storage structure for the back-calculated PCS values. 400 @type pcs_theta: numpy rank-2 array 401 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors. 402 @type pcs_theta_err: numpy rank-2 array 403 @keyword missing_pcs: A structure used to indicate which PCS values are missing. 404 @type missing_pcs: numpy rank-2 array 405 """ 406 407 # Clear the data structures. 408 pcs_theta[:] = 0.0 409 pcs_theta_err[:] = 0.0 410 411 # Fast frame shift. 412 Ri = dot(R_eigen, tensordot(Ri_prime, RT_eigen, axes=1)) 413 Ri = swapaxes(Ri, 0, 1) 414 415 # Unpack the points. 416 theta, phi = points 417 418 # Calculate theta_max. 419 theta_max = tmax_pseudo_ellipse_array(phi, theta_x, theta_y) 420 421 # Loop over the samples. 422 num = 0 423 for i in range(len(points[0])): 424 # The maximum number of points has been reached (well, surpassed by one so exit the loop before it is used). 425 if num == max_points: 426 break 427 428 # As theta_x <= theta_y, check if theta is outside of the isotropic cone defined by theta_y to minimise calculations for speed. 429 if theta[i] > theta_y: 430 continue 431 432 # Outside of the distribution, so skip the point. 433 if theta[i] > theta_max[i]: 434 continue 435 436 # Calculate the PCSs for this state. 437 pcs_pivot_motion_torsionless_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) 438 439 # Increment the number of points. 440 num += 1 441 442 # Default to the rigid state if no points lie in the distribution. 443 if num == 0: 444 # Fast identity frame shift. 445 Ri_prime = eye(3, dtype=float64) 446 Ri = dot(R_eigen, tensordot(Ri_prime, RT_eigen, axes=1)) 447 Ri = swapaxes(Ri, 0, 1) 448 449 # Calculate the PCSs for this state. 450 pcs_pivot_motion_torsionless_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) 451 452 # Multiply the constant. 453 multiply(c, pcs_theta, pcs_theta) 454 455 # Average the PCS. 456 else: 457 multiply(c, pcs_theta, pcs_theta) 458 divide(pcs_theta, float(num), pcs_theta)

459 460

461 -def pcs_numeric_quad_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):

462 """Determine the averaged PCS value via numerical integration. 463 464 @keyword theta_x: The x-axis half cone angle. 465 @type theta_x: float 466 @keyword theta_y: The y-axis half cone angle. 467 @type theta_y: float 468 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units). 469 @type c: float 470 @keyword r_pivot_atom: The pivot point to atom vector. 471 @type r_pivot_atom: numpy rank-1, 3D array 472 @keyword r_ln_pivot: The lanthanide position to pivot point vector. 473 @type r_ln_pivot: numpy rank-1, 3D array 474 @keyword A: The full alignment tensor of the non-moving domain. 475 @type A: numpy rank-2, 3D array 476 @keyword R_eigen: The eigenframe rotation matrix. 477 @type R_eigen: numpy rank-2, 3D array 478 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations). 479 @type RT_eigen: numpy rank-2, 3D array 480 @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. 481 @type Ri_prime: numpy rank-2, 3D array 482 @return: The averaged PCS value. 483 @rtype: float 484 """ 485 486 def pseudo_ellipse(phi): 487 """The pseudo-ellipse wrapper formula.""" 488 489 return tmax_pseudo_ellipse(phi, theta_x, theta_y)

490 491 # Perform numerical integration. 492 result = dblquad(pcs_pivot_motion_torsionless_quad_int, -pi, pi, lambda phi: 0.0, pseudo_ellipse, args=(r_pivot_atom, r_ln_pivot, A, R_eigen, RT_eigen, Ri_prime)) 493 494 # The surface area normalisation factor. 495 SA = pec(theta_x, theta_y) 496 497 # Return the value. 498 return c * result[0] / SA 499

Source Code for Module lib.frame_order.pseudo_ellipse_torsionless