lib.frame_order.double

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 double rotor frame order model.""" 24 25 # Python module imports. 26 from math import cos, pi, sin 27 from numpy import add, divide, dot, eye, float64, multiply, sinc, swapaxes, tensordot 28 try: 29 from scipy.integrate import dblquad 30 except ImportError: 31 pass 32 33 # relax module imports. 34 from lib.compat import norm 35 from lib.frame_order.matrix_ops import rotate_daeg 36 37

38 -def compile_1st_matrix_double_rotor(matrix, R_eigen, smax1, smax2):

39 """Generate the 1st degree Frame Order matrix for the double rotor model. 40 41 @param matrix: The Frame Order matrix, 1st degree to be populated. 42 @type matrix: numpy 3D, rank-2 array 43 @param R_eigen: The eigenframe rotation matrix. 44 @type R_eigen: numpy 3D, rank-2 array 45 @param smax1: The maximum torsion angle for the first rotor. 46 @type smax1: float 47 @param smax2: The maximum torsion angle for the second rotor. 48 @type smax2: float 49 """ 50 51 # Repetitive trig calculations. 52 sinc_smax1 = sinc(smax1/pi) 53 sinc_smax2 = sinc(smax2/pi) 54 55 # Numerical integration of phi of each element. 56 matrix[0, 0] = sinc_smax1 57 matrix[1, 1] = sinc_smax2 58 matrix[2, 2] = sinc_smax1 * sinc_smax2 59 60 # Rotate and return the frame order matrix. 61 return rotate_daeg(matrix, R_eigen)

62 63

64 -def compile_2nd_matrix_double_rotor(matrix, Rx2_eigen, smax1, smax2):

65 """Generate the rotated 2nd degree Frame Order matrix for the double rotor model. 66 67 The cone axis is assumed to be parallel to the z-axis in the eigenframe. 68 69 70 @param matrix: The Frame Order matrix, 2nd degree to be populated. 71 @type matrix: numpy 9D, rank-2 array 72 @param Rx2_eigen: The Kronecker product of the eigenframe rotation matrix with itself. 73 @type Rx2_eigen: numpy 9D, rank-2 array 74 @param smax1: The maximum torsion angle for the first rotor. 75 @type smax1: float 76 @param smax2: The maximum torsion angle for the second rotor. 77 @type smax2: float 78 """ 79 80 # Zeros. 81 matrix[:] = 0.0 82 83 # Repetitive trig calculations. 84 sinc_smax1 = sinc(smax1/pi) 85 sinc_2smax1 = sinc(2.0*smax1/pi) 86 sinc_2smax1p1 = sinc_2smax1 + 1.0 87 sinc_2smax1n1 = sinc_2smax1 - 1.0 88 sinc_smax2 = sinc(smax2/pi) 89 sinc_2smax2 = sinc(2.0*smax2/pi) 90 sinc_2smax2p1 = sinc_2smax2 + 1.0 91 sinc_2smax2n1 = sinc_2smax2 - 1.0 92 93 # Diagonal. 94 matrix[0, 0] = sinc_2smax1 + 1.0 95 matrix[1, 1] = 2.0 * sinc_smax1 * sinc_smax2 96 matrix[2, 2] = sinc_smax2 * sinc_2smax1p1 97 matrix[3, 3] = matrix[1, 1] 98 matrix[4, 4] = sinc_2smax2p1 99 matrix[5, 5] = sinc_smax1 * sinc_2smax2p1 100 matrix[6, 6] = matrix[2, 2] 101 matrix[7, 7] = matrix[5, 5] 102 matrix[8, 8] = 0.5 * sinc_2smax1p1 * sinc_2smax2p1 103 104 # Off diagonal set 1. 105 matrix[4, 0] = 0.5 * sinc_2smax1n1 * sinc_2smax2n1 106 matrix[0, 8] = -sinc_2smax1n1 107 matrix[8, 0] = -0.5 * sinc_2smax1n1 * sinc_2smax2p1 108 matrix[4, 8] = -0.5 * sinc_2smax1p1 * sinc_2smax2n1 109 matrix[8, 4] = -sinc_2smax2n1 110 111 # Off diagonal set 2. 112 matrix[2, 6] = matrix[6, 2] = sinc_smax2 * sinc_2smax1n1 113 matrix[5, 7] = matrix[7, 5] = sinc_smax1 * sinc_2smax2n1 114 115 # Divide by 2. 116 multiply(0.5, matrix, matrix) 117 118 # Rotate and return the frame order matrix. 119 return rotate_daeg(matrix, Rx2_eigen)

120 121

122 -def pcs_numeric_qr_int_double_rotor(points=None, max_points=None, sigma_max=None, sigma_max_2=None, c=None, full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None, r_ln_pivot=None, r_inter_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None, Ri2_prime=None, pcs_theta=None, pcs_theta_err=None, missing_pcs=None):

123 """The averaged PCS value via numerical integration for the double rotor frame order model. 124 125 @keyword points: The Sobol points in the torsion-tilt angle space. 126 @type points: numpy rank-2, 3D array 127 @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. 128 @type max_points: int 129 @keyword sigma_max: The maximum opening angle for the first rotor. 130 @type sigma_max: float 131 @keyword sigma_max_2: The maximum opening angle for the second rotor. 132 @type sigma_max_2: float 133 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units). 134 @type c: numpy rank-1 array 135 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 136 @type full_in_ref_frame: numpy rank-1 array 137 @keyword r_pivot_atom: The pivot point to atom vector. 138 @type r_pivot_atom: numpy rank-2, 3D array 139 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector. 140 @type r_pivot_atom_rev: numpy rank-2, 3D array 141 @keyword r_ln_pivot: The lanthanide position to pivot point vector. 142 @type r_ln_pivot: numpy rank-2, 3D array 143 @keyword r_inter_pivot: The vector between the two pivots. 144 @type r_inter_pivot: numpy rank-1, 3D array 145 @keyword A: The full alignment tensor of the non-moving domain. 146 @type A: numpy rank-2, 3D array 147 @keyword R_eigen: The eigenframe rotation matrix. 148 @type R_eigen: numpy rank-2, 3D array 149 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations). 150 @type RT_eigen: numpy rank-2, 3D array 151 @keyword Ri_prime: The array of pre-calculated rotation matrices for the in-frame double rotor motion for the 1st mode of motion, used to calculate the PCS for each state i in the numerical integration. 152 @type Ri_prime: numpy rank-3, array of 3D arrays 153 @keyword Ri2_prime: The array of pre-calculated rotation matrices for the in-frame double rotor motion for the 2nd mode of motion, used to calculate the PCS for each state i in the numerical integration. 154 @type Ri2_prime: numpy rank-3, array of 3D arrays 155 @keyword pcs_theta: The storage structure for the back-calculated PCS values. 156 @type pcs_theta: numpy rank-2 array 157 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors. 158 @type pcs_theta_err: numpy rank-2 array 159 @keyword missing_pcs: A structure used to indicate which PCS values are missing. 160 @type missing_pcs: numpy rank-2 array 161 """ 162 163 # Clear the data structures. 164 pcs_theta[:] = 0.0 165 pcs_theta_err[:] = 0.0 166 167 # Fast frame shift. 168 Ri = dot(R_eigen, tensordot(Ri_prime, RT_eigen, axes=1)) 169 Ri = swapaxes(Ri, 0, 1) 170 Ri2 = dot(R_eigen, tensordot(Ri2_prime, RT_eigen, axes=1)) 171 Ri2 = swapaxes(Ri2, 0, 1) 172 173 # Unpack the points. 174 sigma, sigma2 = points 175 176 # Loop over the samples. 177 num = 0 178 for i in range(len(points[0])): 179 # The maximum number of points has been reached (well, surpassed by one so exit the loop before it is used). 180 if num == max_points: 181 break 182 183 # Outside of the distribution, so skip the point. 184 if abs(sigma[i]) > sigma_max: 185 continue 186 if abs(sigma2[i]) > sigma_max_2: 187 continue 188 189 # Calculate the PCSs for this state. 190 pcs_pivot_motion_double_rotor_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, r_inter_pivot=r_inter_pivot, A=A, Ri=Ri[i], Ri2=Ri2[i], pcs_theta=pcs_theta, pcs_theta_err=pcs_theta_err, missing_pcs=missing_pcs) 191 192 # Increment the number of points. 193 num += 1 194 195 # Default to the rigid state if no points lie in the distribution. 196 if num == 0: 197 # Fast identity frame shift. 198 Ri_prime = eye(3, dtype=float64) 199 Ri = dot(R_eigen, tensordot(Ri_prime, RT_eigen, axes=1)) 200 Ri = swapaxes(Ri, 0, 1) 201 202 # Calculate the PCSs for this state. 203 pcs_pivot_motion_double_rotor_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, r_inter_pivot=r_inter_pivot, A=A, Ri=Ri, Ri2=Ri, pcs_theta=pcs_theta, pcs_theta_err=pcs_theta_err, missing_pcs=missing_pcs) 204 205 # Multiply the constant. 206 multiply(c, pcs_theta, pcs_theta) 207 208 # Average the PCS. 209 else: 210 multiply(c, pcs_theta, pcs_theta) 211 divide(pcs_theta, float(num), pcs_theta)

212 213

214 -def pcs_numeric_quad_int_double_rotor(sigma_max=None, sigma_max_2=None, c=None, r_pivot_atom=None, r_ln_pivot=None, r_inter_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None, Ri2_prime=None):

215 """Determine the averaged PCS value via SciPy quadratic numerical integration. 216 217 @keyword sigma_max: The maximum opening angle for the first rotor. 218 @type sigma_max: float 219 @keyword sigma_max_2: The maximum opening angle for the second rotor. 220 @type sigma_max_2: float 221 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units). 222 @type c: numpy rank-1 array 223 @keyword r_pivot_atom: The pivot point to atom vector. 224 @type r_pivot_atom: numpy rank-1, 3D array 225 @keyword r_ln_pivot: The lanthanide position to pivot point vector. 226 @type r_ln_pivot: numpy rank-1, 3D array 227 @keyword r_inter_pivot: The vector between the two pivots. 228 @type r_inter_pivot: numpy rank-1, 3D array 229 @keyword A: The full alignment tensor of the non-moving domain. 230 @type A: numpy rank-2, 3D array 231 @keyword R_eigen: The eigenframe rotation matrix. 232 @type R_eigen: numpy rank-2, 3D array 233 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations). 234 @type RT_eigen: numpy rank-2, 3D array 235 @keyword Ri_prime: The array of pre-calculated rotation matrices for the in-frame double rotor motion for the 1st mode of motion, used to calculate the PCS for each state i in the numerical integration. 236 @type Ri_prime: numpy rank-2, 3D array 237 @keyword Ri2_prime: The array of pre-calculated rotation matrices for the in-frame double rotor motion for the 2nd mode of motion, used to calculate the PCS for each state i in the numerical integration. 238 @type Ri2_prime: numpy rank-2, 3D array 239 """ 240 241 # Preset the 1st rotation matrix elements for state i. 242 Ri_prime[0, 1] = 0.0 243 Ri_prime[1, 0] = 0.0 244 Ri_prime[1, 1] = 1.0 245 Ri_prime[1, 2] = 0.0 246 Ri_prime[2, 1] = 0.0 247 248 # Preset the 2nd rotation matrix elements for state i. 249 Ri2_prime[0, 0] = 1.0 250 Ri2_prime[0, 1] = 0.0 251 Ri2_prime[0, 2] = 0.0 252 Ri2_prime[1, 0] = 0.0 253 Ri2_prime[2, 0] = 0.0 254 255 # Perform numerical integration. 256 result = dblquad(pcs_pivot_motion_double_rotor_quad_int, -sigma_max, sigma_max, lambda sigma2: -sigma_max_2, lambda sigma2: sigma_max_2, args=(r_pivot_atom, r_ln_pivot, r_inter_pivot, A, R_eigen, RT_eigen, Ri_prime, Ri2_prime)) 257 258 # The surface area normalisation factor. 259 SA = 4.0 * sigma_max * sigma_max_2 260 261 # Return the value. 262 return c * result[0] / SA

263 264

265 -def pcs_pivot_motion_double_rotor_qr_int(full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None, r_ln_pivot=None, r_inter_pivot=None, A=None, Ri=None, Ri2=None, pcs_theta=None, pcs_theta_err=None, missing_pcs=None):

266 """Calculate the PCS value after a pivoted motion for the double rotor model. 267 268 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 269 @type full_in_ref_frame: numpy rank-1 array 270 @keyword r_pivot_atom: The pivot point to atom vector. 271 @type r_pivot_atom: numpy rank-2, 3D array 272 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector. 273 @type r_pivot_atom_rev: numpy rank-2, 3D array 274 @keyword r_ln_pivot: The lanthanide position to pivot point vector. 275 @type r_ln_pivot: numpy rank-2, 3D array 276 @keyword r_inter_pivot: The vector between the two pivots. 277 @type r_inter_pivot: numpy rank-1, 3D array 278 @keyword A: The full alignment tensor of the non-moving domain. 279 @type A: numpy rank-2, 3D array 280 @keyword Ri: The frame-shifted, pre-calculated rotation matrix for state i for the 1st mode of motion. 281 @type Ri: numpy rank-2, 3D array 282 @keyword Ri2: The frame-shifted, pre-calculated rotation matrix for state i for the 2nd mode of motion. 283 @type Ri2: numpy rank-2, 3D array 284 @keyword pcs_theta: The storage structure for the back-calculated PCS values. 285 @type pcs_theta: numpy rank-2 array 286 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors. 287 @type pcs_theta_err: numpy rank-2 array 288 @keyword missing_pcs: A structure used to indicate which PCS values are missing. 289 @type missing_pcs: numpy rank-2 array 290 """ 291 292 # Rotate the first pivot to atomic position vectors. 293 rot_vect = dot(r_pivot_atom, Ri) 294 295 # Add the inter-pivot vector to obtain the 2nd pivot to atomic position vectors. 296 add(r_inter_pivot, rot_vect, rot_vect) 297 298 # Rotate the 2nd pivot to atomic position vectors. 299 rot_vect = dot(rot_vect, Ri2) 300 301 # Add the lanthanide to pivot vector. 302 add(rot_vect, r_ln_pivot, rot_vect) 303 304 # The vector length (to the 5th power). 305 length = 1.0 / norm(rot_vect, axis=1)**5 306 307 # The reverse vectors and lengths. 308 if min(full_in_ref_frame) == 0: 309 rot_vect_rev = dot(r_pivot_atom_rev, Ri) 310 add(r_inter_pivot, rot_vect_rev, rot_vect_rev) 311 rot_vect_rev = dot(rot_vect_rev, Ri2) 312 add(rot_vect_rev, r_ln_pivot, rot_vect_rev) 313 length_rev = 1.0 / norm(rot_vect_rev, axis=1)**5 314 315 # Loop over the atoms. 316 for j in range(len(r_pivot_atom[:, 0])): 317 # Loop over the alignments. 318 for i in range(len(pcs_theta)): 319 # Skip missing data. 320 if missing_pcs[i, j]: 321 continue 322 323 # The projection. 324 if full_in_ref_frame[i]: 325 proj = dot(rot_vect[j], dot(A[i], rot_vect[j])) 326 length_i = length[j] 327 else: 328 proj = dot(rot_vect_rev[j], dot(A[i], rot_vect_rev[j])) 329 length_i = length_rev[j] 330 331 # The PCS. 332 pcs_theta[i, j] += proj * length_i

333 334

335 -def pcs_pivot_motion_double_rotor_quad_int(sigma_i, sigma2_i, r_pivot_atom, r_ln_pivot, r_inter_pivot, A, R_eigen, RT_eigen, Ri_prime, Ri2_prime):

336 """Calculate the PCS value after a pivoted motion for the double rotor model. 337 338 @param sigma_i: The 1st torsion angle for state i. 339 @type sigma_i: float 340 @param sigma2_i: The 1st torsion angle for state i. 341 @type sigma2_i: float 342 @param r_pivot_atom: The pivot point to atom vector. 343 @type r_pivot_atom: numpy rank-2, 3D array 344 @param r_ln_pivot: The lanthanide position to pivot point vector. 345 @type r_ln_pivot: numpy rank-2, 3D array 346 @param r_inter_pivot: The vector between the two pivots. 347 @type r_inter_pivot: numpy rank-1, 3D array 348 @param A: The full alignment tensor of the non-moving domain. 349 @type A: numpy rank-2, 3D array 350 @param R_eigen: The eigenframe rotation matrix. 351 @type R_eigen: numpy rank-2, 3D array 352 @param RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations). 353 @type RT_eigen: numpy rank-2, 3D array 354 @param Ri_prime: The empty rotation matrix for state i. 355 @type Ri_prime: numpy rank-2, 3D array 356 @param Ri2_prime: The 2nd empty rotation matrix for state i. 357 @type Ri2_prime: numpy rank-2, 3D array 358 @return: The PCS value for the changed position. 359 @rtype: float 360 """ 361 362 # The 1st rotation matrix. 363 c_sigma = cos(sigma_i) 364 s_sigma = sin(sigma_i) 365 Ri_prime[0, 0] = c_sigma 366 Ri_prime[0, 2] = s_sigma 367 Ri_prime[2, 0] = -s_sigma 368 Ri_prime[2, 2] = c_sigma 369 370 # The 2nd rotation matrix. 371 c_sigma = cos(sigma2_i) 372 s_sigma = sin(sigma2_i) 373 Ri2_prime[1, 1] = c_sigma 374 Ri2_prime[1, 2] = -s_sigma 375 Ri2_prime[2, 1] = s_sigma 376 Ri2_prime[2, 2] = c_sigma 377 378 # The rotations. 379 Ri = dot(R_eigen, dot(Ri_prime, RT_eigen)) 380 Ri2 = dot(R_eigen, dot(Ri2_prime, RT_eigen)) 381 382 # Rotate the first pivot to atomic position vectors. 383 rot_vect = dot(r_pivot_atom, Ri) 384 385 # Add the inter-pivot vector to obtain the 2nd pivot to atomic position vectors. 386 add(r_inter_pivot, rot_vect, rot_vect) 387 388 # Rotate the 2nd pivot to atomic position vectors. 389 rot_vect = dot(rot_vect, Ri2) 390 391 # Add the lanthanide to pivot vector. 392 add(rot_vect, r_ln_pivot, rot_vect) 393 394 # The vector length. 395 length = norm(rot_vect) 396 397 # The projection. 398 proj = dot(rot_vect, dot(A, rot_vect)) 399 400 # The PCS. 401 pcs = proj / length**5 402 403 # Return the PCS value (without the PCS constant). 404 return pcs

405

Source Code for Module lib.frame_order.double_rotor