lib.frame_order.pseudo_ellipse

1 ############################################################################### 2 # # 3 # Copyright (C) 2009-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, sqrt 27 try: 28 from scipy.integrate import dblquad, quad 29 except ImportError: 30 pass 31 32 # relax module imports. 33 from lib.geometry.pec import pec 34 from lib.frame_order.matrix_ops import pcs_pivot_motion_torsionless, pcs_pivot_motion_torsionless_qrint, rotate_daeg 35 from lib.frame_order.pseudo_ellipse import tmax_pseudo_ellipse 36 37

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

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

77 78

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

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

97 98

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

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

117 118

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

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

137 138

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

140 """The theta partial integral of the 2nd degree Frame Order matrix for the torsionless pseudo-ellipse. 141 142 @param phi: The azimuthal tilt-torsion angle. 143 @type phi: float 144 @param x: The cone opening angle along x. 145 @type x: float 146 @param y: The cone opening angle along y. 147 @type y: float 148 @return: The theta partial integral. 149 @rtype: float 150 """ 151 152 # Theta max. 153 tmax = tmax_pseudo_ellipse(phi, x, y) 154 155 # The theta integral. 156 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)

157 158

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

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

177 178

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

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

197 198

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

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

217 218

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

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

237 238

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

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

257 258

259 -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):

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

288 289 # Perform numerical integration. 290 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)) 291 292 # The surface area normalisation factor. 293 SA = pec(theta_x, theta_y) 294 295 # Return the value. 296 return c * result[0] / SA 297 298

299 -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):

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

372

Source Code for Module lib.frame_order.pseudo_ellipse_torsionless