lib.frame_order.matrix

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 # Python module imports. 26 from math import cos, sin 27 from numpy import dot, transpose 28 from numpy.linalg import norm 29 30 # relax module imports. 31 from lib.compat import norm 32 from lib.linear_algebra.kronecker_product import transpose_23 33 34

35 -def daeg_to_rotational_superoperator(daeg, Rsuper):

36 """Convert the frame order matrix (daeg) to the rotational superoperator. 37 38 @param daeg: The second degree frame order matrix, daeg. This must be in the Kronecker product layout. 39 @type daeg: numpy 9D, rank-2 array or numpy 3D, rank-4 array 40 @param Rsuper: The rotational superoperator structure to be populated. 41 @type Rsuper: numpy 5D, rank-2 array 42 """ 43 44 # First perform the T23 transpose. 45 transpose_23(daeg) 46 47 # Convert to rank-4. 48 orig_shape = daeg.shape 49 daeg.shape = (3, 3, 3, 3) 50 51 # First column of the superoperator. 52 Rsuper[0, 0] = daeg[0, 0, 0, 0] - daeg[2, 0, 2, 0] 53 Rsuper[1, 0] = daeg[0, 1, 0, 1] - daeg[2, 1, 2, 1] 54 Rsuper[2, 0] = daeg[0, 0, 0, 1] - daeg[2, 0, 2, 1] 55 Rsuper[3, 0] = daeg[0, 0, 0, 2] - daeg[2, 0, 2, 2] 56 Rsuper[4, 0] = daeg[0, 1, 0, 2] - daeg[2, 1, 2, 2] 57 58 # Second column of the superoperator. 59 Rsuper[0, 1] = daeg[1, 0, 1, 0] - daeg[2, 0, 2, 0] 60 Rsuper[1, 1] = daeg[1, 1, 1, 1] - daeg[2, 1, 2, 1] 61 Rsuper[2, 1] = daeg[1, 0, 1, 1] - daeg[2, 0, 2, 1] 62 Rsuper[3, 1] = daeg[1, 0, 1, 2] - daeg[2, 0, 2, 2] 63 Rsuper[4, 1] = daeg[1, 1, 1, 2] - daeg[2, 1, 2, 2] 64 65 # Third column of the superoperator. 66 Rsuper[0, 2] = daeg[0, 0, 1, 0] + daeg[1, 0, 0, 0] 67 Rsuper[1, 2] = daeg[0, 1, 1, 1] + daeg[1, 1, 0, 1] 68 Rsuper[2, 2] = daeg[0, 0, 1, 1] + daeg[1, 0, 0, 1] 69 Rsuper[3, 2] = daeg[0, 0, 1, 2] + daeg[1, 0, 0, 2] 70 Rsuper[4, 2] = daeg[0, 1, 1, 2] + daeg[1, 1, 0, 2] 71 72 # Fourth column of the superoperator. 73 Rsuper[0, 3] = daeg[0, 0, 2, 0] + daeg[2, 0, 0, 0] 74 Rsuper[1, 3] = daeg[0, 1, 2, 1] + daeg[2, 1, 0, 1] 75 Rsuper[2, 3] = daeg[0, 0, 2, 1] + daeg[2, 0, 0, 1] 76 Rsuper[3, 3] = daeg[0, 0, 2, 2] + daeg[2, 0, 0, 2] 77 Rsuper[4, 3] = daeg[0, 1, 2, 2] + daeg[2, 1, 0, 2] 78 79 # Fifth column of the superoperator. 80 Rsuper[0, 4] = daeg[1, 0, 2, 0] + daeg[2, 0, 1, 0] 81 Rsuper[1, 4] = daeg[1, 1, 2, 1] + daeg[2, 1, 1, 1] 82 Rsuper[2, 4] = daeg[1, 0, 2, 1] + daeg[2, 0, 1, 1] 83 Rsuper[3, 4] = daeg[1, 0, 2, 2] + daeg[2, 0, 1, 2] 84 Rsuper[4, 4] = daeg[1, 1, 2, 2] + daeg[2, 1, 1, 2] 85 86 # Revert the shape. 87 daeg.shape = orig_shape 88 89 # Undo the T23 transpose. 90 transpose_23(daeg)

91 92

93 -def pcs_pivot_motion_full_qr_int(full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None, r_ln_pivot=None, A=None, Ri=None, pcs_theta=None, pcs_theta_err=None, missing_pcs=None):

94 """Calculate the PCS value after a pivoted motion for the isotropic cone model. 95 96 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 97 @type full_in_ref_frame: numpy rank-1 array 98 @keyword r_pivot_atom: The pivot point to atom vector. 99 @type r_pivot_atom: numpy rank-2, 3D array 100 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector. 101 @type r_pivot_atom_rev: numpy rank-2, 3D array 102 @keyword r_ln_pivot: The lanthanide position to pivot point vector. 103 @type r_ln_pivot: numpy rank-2, 3D array 104 @keyword A: The full alignment tensor of the non-moving domain. 105 @type A: numpy rank-2, 3D array 106 @keyword Ri: The frame-shifted, pre-calculated rotation matrix for state i. 107 @type Ri: numpy rank-2, 3D array 108 @keyword pcs_theta: The storage structure for the back-calculated PCS values. 109 @type pcs_theta: numpy rank-2 array 110 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors. 111 @type pcs_theta_err: numpy rank-2 array 112 @keyword missing_pcs: A structure used to indicate which PCS values are missing. 113 @type missing_pcs: numpy rank-2 array 114 """ 115 116 # Pre-calculate all the new vectors. 117 rot_vect = dot(r_pivot_atom, Ri) + r_ln_pivot 118 119 # The vector length (to the 5th power). 120 length = 1.0 / norm(rot_vect, axis=1)**5 121 122 # The reverse vectors and lengths. 123 if min(full_in_ref_frame) == 0: 124 rot_vect_rev = dot(r_pivot_atom_rev, Ri) + r_ln_pivot 125 length_rev = 1.0 / norm(rot_vect_rev, axis=1)**5 126 127 # Loop over the atoms. 128 for j in range(len(r_pivot_atom[:, 0])): 129 # Loop over the alignments. 130 for i in range(len(pcs_theta)): 131 # Skip missing data. 132 if missing_pcs[i, j]: 133 continue 134 135 # The projection. 136 if full_in_ref_frame[i]: 137 proj = dot(rot_vect[j], dot(A[i], rot_vect[j])) 138 length_i = length[j] 139 else: 140 proj = dot(rot_vect_rev[j], dot(A[i], rot_vect_rev[j])) 141 length_i = length_rev[j] 142 143 # The PCS. 144 pcs_theta[i, j] += proj * length_i

145 146

147 -def pcs_pivot_motion_full_quad_int(theta_i, phi_i, sigma_i, r_pivot_atom, r_ln_pivot, A, R_eigen, RT_eigen, Ri_prime):

148 """Calculate the PCS value after a pivoted motion for the isotropic cone model. 149 150 @param theta_i: The half cone opening angle (polar angle). 151 @type theta_i: float 152 @param phi_i: The cone azimuthal angle. 153 @type phi_i: float 154 @param sigma_i: The torsion angle for state i. 155 @type sigma_i: float 156 @param r_pivot_atom: The pivot point to atom vector. 157 @type r_pivot_atom: numpy rank-1, 3D array 158 @param r_ln_pivot: The lanthanide position to pivot point vector. 159 @type r_ln_pivot: numpy rank-1, 3D array 160 @param A: The full alignment tensor of the non-moving domain. 161 @type A: numpy rank-2, 3D array 162 @param R_eigen: The eigenframe rotation matrix. 163 @type R_eigen: numpy rank-2, 3D array 164 @param RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations). 165 @type RT_eigen: numpy rank-2, 3D array 166 @param Ri_prime: The empty rotation matrix for the in-frame isotropic cone motion for state i. 167 @type Ri_prime: numpy rank-2, 3D array 168 @return: The PCS value for the changed position. 169 @rtype: float 170 """ 171 172 # The rotation matrix. 173 c_theta = cos(theta_i) 174 s_theta = sin(theta_i) 175 c_phi = cos(phi_i) 176 s_phi = sin(phi_i) 177 c_sigma_phi = cos(sigma_i - phi_i) 178 s_sigma_phi = sin(sigma_i - phi_i) 179 c_phi_c_theta = c_phi * c_theta 180 s_phi_c_theta = s_phi * c_theta 181 Ri_prime[0, 0] = c_phi_c_theta*c_sigma_phi - s_phi*s_sigma_phi 182 Ri_prime[0, 1] = -c_phi_c_theta*s_sigma_phi - s_phi*c_sigma_phi 183 Ri_prime[0, 2] = c_phi*s_theta 184 Ri_prime[1, 0] = s_phi_c_theta*c_sigma_phi + c_phi*s_sigma_phi 185 Ri_prime[1, 1] = -s_phi_c_theta*s_sigma_phi + c_phi*c_sigma_phi 186 Ri_prime[1, 2] = s_phi*s_theta 187 Ri_prime[2, 0] = -s_theta*c_sigma_phi 188 Ri_prime[2, 1] = s_theta*s_sigma_phi 189 Ri_prime[2, 2] = c_theta 190 191 # The rotation. 192 R_i = dot(R_eigen, dot(Ri_prime, RT_eigen)) 193 194 # Calculate the new vector. 195 vect = dot(R_i, r_pivot_atom) + r_ln_pivot 196 197 # The vector length. 198 length = norm(vect) 199 200 # The projection. 201 proj = dot(vect, dot(A, vect)) 202 203 # The PCS (with sine surface normalisation). 204 pcs = proj / length**5 * s_theta 205 206 # Return the PCS value (without the PCS constant). 207 return pcs

208 209

210 -def pcs_pivot_motion_torsionless_qr_int(full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None, r_ln_pivot=None, A=None, Ri=None, pcs_theta=None, pcs_theta_err=None, missing_pcs=None):

211 """Calculate the PCS value after a pivoted motion for the isotropic cone model. 212 213 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 214 @type full_in_ref_frame: numpy rank-1 array 215 @keyword r_pivot_atom: The pivot point to atom vector. 216 @type r_pivot_atom: numpy rank-2, 3D array 217 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector. 218 @type r_pivot_atom_rev: numpy rank-2, 3D array 219 @keyword r_ln_pivot: The lanthanide position to pivot point vector. 220 @type r_ln_pivot: numpy rank-2, 3D array 221 @keyword A: The full alignment tensor of the non-moving domain. 222 @type A: numpy rank-2, 3D array 223 @keyword Ri: The frame-shifted, pre-calculated rotation matrix for state i. 224 @type Ri: numpy rank-2, 3D array 225 @keyword pcs_theta: The storage structure for the back-calculated PCS values. 226 @type pcs_theta: numpy rank-2 array 227 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors. 228 @type pcs_theta_err: numpy rank-2 array 229 @keyword missing_pcs: A structure used to indicate which PCS values are missing. 230 @type missing_pcs: numpy rank-2 array 231 """ 232 233 # Pre-calculate all the new vectors. 234 rot_vect = dot(r_pivot_atom, Ri) + r_ln_pivot 235 236 # The vector length (to the 5th power). 237 length = 1.0 / norm(rot_vect, axis=1)**5 238 239 # The reverse vectors and lengths. 240 if min(full_in_ref_frame) == 0: 241 rot_vect_rev = dot(r_pivot_atom_rev, Ri) + r_ln_pivot 242 length_rev = 1.0 / norm(rot_vect_rev, axis=1)**5 243 244 # Loop over the atoms. 245 for j in range(len(r_pivot_atom[:, 0])): 246 # Loop over the alignments. 247 for i in range(len(pcs_theta)): 248 # Skip missing data. 249 if missing_pcs[i, j]: 250 continue 251 252 # The projection. 253 if full_in_ref_frame[i]: 254 proj = dot(rot_vect[j], dot(A[i], rot_vect[j])) 255 length_i = length[j] 256 else: 257 proj = dot(rot_vect_rev[j], dot(A[i], rot_vect_rev[j])) 258 length_i = length_rev[j] 259 260 # The PCS. 261 pcs_theta[i, j] += proj * length_i

262 263

264 -def pcs_pivot_motion_torsionless_quad_int(theta_i, phi_i, r_pivot_atom, r_ln_pivot, A, R_eigen, RT_eigen, Ri_prime):

265 """Calculate the PCS value after a pivoted motion for the isotropic cone model. 266 267 @param theta_i: The half cone opening angle (polar angle). 268 @type theta_i: float 269 @param phi_i: The cone azimuthal angle. 270 @type phi_i: float 271 @param r_pivot_atom: The pivot point to atom vector. 272 @type r_pivot_atom: numpy rank-1, 3D array 273 @param r_ln_pivot: The lanthanide position to pivot point vector. 274 @type r_ln_pivot: numpy rank-1, 3D array 275 @param A: The full alignment tensor of the non-moving domain. 276 @type A: numpy rank-2, 3D array 277 @param R_eigen: The eigenframe rotation matrix. 278 @type R_eigen: numpy rank-2, 3D array 279 @param RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations). 280 @type RT_eigen: numpy rank-2, 3D array 281 @param Ri_prime: The empty rotation matrix for the in-frame isotropic cone motion for state i. 282 @type Ri_prime: numpy rank-2, 3D array 283 @return: The PCS value for the changed position. 284 @rtype: float 285 """ 286 287 # The rotation matrix. 288 c_theta = cos(theta_i) 289 s_theta = sin(theta_i) 290 c_phi = cos(phi_i) 291 s_phi = sin(phi_i) 292 c_phi_c_theta = c_phi * c_theta 293 s_phi_c_theta = s_phi * c_theta 294 Ri_prime[0, 0] = c_phi_c_theta*c_phi + s_phi**2 295 Ri_prime[0, 1] = c_phi_c_theta*s_phi - c_phi*s_phi 296 Ri_prime[0, 2] = c_phi*s_theta 297 Ri_prime[1, 0] = s_phi_c_theta*c_phi - c_phi*s_phi 298 Ri_prime[1, 1] = s_phi_c_theta*s_phi + c_phi**2 299 Ri_prime[1, 2] = s_phi*s_theta 300 Ri_prime[2, 0] = -s_theta*c_phi 301 Ri_prime[2, 1] = -s_theta*s_phi 302 Ri_prime[2, 2] = c_theta 303 304 # The rotation. 305 R_i = dot(R_eigen, dot(Ri_prime, RT_eigen)) 306 307 # Calculate the new vector. 308 vect = dot(R_i, r_pivot_atom) + r_ln_pivot 309 310 # The vector length. 311 length = norm(vect) 312 313 # The projection. 314 proj = dot(vect, dot(A, vect)) 315 316 # The PCS (with sine surface normalisation). 317 pcs = proj / length**5 * s_theta 318 319 # Return the PCS value (without the PCS constant). 320 return pcs

321 322

323 -def reduce_alignment_tensor(D, A, red_tensor):

324 """Calculate the reduction in the alignment tensor caused by the Frame Order matrix. 325 326 This is both the forward rotation notation and Kronecker product arrangement. 327 328 @param D: The Frame Order matrix, 2nd degree to be populated. 329 @type D: numpy 9D, rank-2 array 330 @param A: The full alignment tensor in {Axx, Ayy, Axy, Axz, Ayz} notation. 331 @type A: numpy 5D, rank-1 array 332 @param red_tensor: The structure in {Axx, Ayy, Axy, Axz, Ayz} notation to place the reduced 333 alignment tensor. 334 @type red_tensor: numpy 5D, rank-1 array 335 """ 336 337 # The reduced tensor element A0. 338 red_tensor[0] = (D[0, 0] - D[8, 0])*A[0] 339 red_tensor[0] = red_tensor[0] + (D[4, 0] - D[8, 0])*A[1] 340 red_tensor[0] = red_tensor[0] + (D[1, 0] + D[3, 0])*A[2] 341 red_tensor[0] = red_tensor[0] + (D[2, 0] + D[6, 0])*A[3] 342 red_tensor[0] = red_tensor[0] + (D[5, 0] + D[7, 0])*A[4] 343 344 # The reduced tensor element A1. 345 red_tensor[1] = (D[0, 4] - D[8, 4])*A[0] 346 red_tensor[1] = red_tensor[1] + (D[4, 4] - D[8, 4])*A[1] 347 red_tensor[1] = red_tensor[1] + (D[1, 4] + D[3, 4])*A[2] 348 red_tensor[1] = red_tensor[1] + (D[2, 4] + D[6, 4])*A[3] 349 red_tensor[1] = red_tensor[1] + (D[5, 4] + D[7, 4])*A[4] 350 351 # The reduced tensor element A2. 352 red_tensor[2] = (D[0, 1] - D[8, 1])*A[0] 353 red_tensor[2] = red_tensor[2] + (D[4, 1] - D[8, 1])*A[1] 354 red_tensor[2] = red_tensor[2] + (D[1, 1] + D[3, 1])*A[2] 355 red_tensor[2] = red_tensor[2] + (D[2, 1] + D[6, 1])*A[3] 356 red_tensor[2] = red_tensor[2] + (D[5, 1] + D[7, 1])*A[4] 357 358 # The reduced tensor element A3. 359 red_tensor[3] = (D[0, 2] - D[8, 2])*A[0] 360 red_tensor[3] = red_tensor[3] + (D[4, 2] - D[8, 2])*A[1] 361 red_tensor[3] = red_tensor[3] + (D[1, 2] + D[3, 2])*A[2] 362 red_tensor[3] = red_tensor[3] + (D[2, 2] + D[6, 2])*A[3] 363 red_tensor[3] = red_tensor[3] + (D[5, 2] + D[7, 2])*A[4] 364 365 # The reduced tensor element A4. 366 red_tensor[4] = (D[0, 5] - D[8, 5])*A[0] 367 red_tensor[4] = red_tensor[4] + (D[4, 5] - D[8, 5])*A[1] 368 red_tensor[4] = red_tensor[4] + (D[1, 5] + D[3, 5])*A[2] 369 red_tensor[4] = red_tensor[4] + (D[2, 5] + D[6, 5])*A[3] 370 red_tensor[4] = red_tensor[4] + (D[5, 5] + D[7, 5])*A[4]

371 372

373 -def reduce_alignment_tensor_symmetric(D, A, red_tensor):

374 """Calculate the reduction in the alignment tensor caused by the Frame Order matrix. 375 376 This is both the forward rotation notation and Kronecker product arrangement. This simplification is due to the symmetry in motion of the pseudo-elliptic and isotropic cones. All element of the frame order matrix where an index appears only once are zero. 377 378 @param D: The Frame Order matrix, 2nd degree to be populated. 379 @type D: numpy 9D, rank-2 array 380 @param A: The full alignment tensor in {Axx, Ayy, Axy, Axz, Ayz} notation. 381 @type A: numpy 5D, rank-1 array 382 @param red_tensor: The structure in {Axx, Ayy, Axy, Axz, Ayz} notation to place the reduced 383 alignment tensor. 384 @type red_tensor: numpy 5D, rank-1 array 385 """ 386 387 # The reduced tensor elements. 388 red_tensor[0] = (D[0, 0] - D[8, 0])*A[0] + (D[4, 0] - D[8, 0])*A[1] 389 red_tensor[1] = (D[0, 4] - D[8, 4])*A[0] + (D[4, 4] - D[8, 4])*A[1] 390 red_tensor[2] = (D[1, 1] + D[3, 1])*A[2] 391 red_tensor[3] = (D[2, 2] + D[6, 2])*A[3] 392 red_tensor[4] = (D[5, 5] + D[7, 5])*A[4]

393 394

395 -def rotate_daeg(matrix, R_eigen):

396 """Rotate the given frame order matrix. 397 398 It is assumed that the frame order matrix is in the Kronecker product form. 399 400 401 @param matrix: The Frame Order matrix to be populated. All degrees are handled. 402 @type matrix: numpy rank-2 array 403 @param R_eigen: The Kronecker product of the eigenframe rotation matrix with itself. 404 @type R_eigen: numpy rank-2 array 405 """ 406 407 # Rotate. 408 matrix_rot = dot(R_eigen, dot(matrix, transpose(R_eigen))) 409 410 # Return the matrix. 411 return matrix_rot

412 413

414 -class Data:

415 """A data container stored in the memo objects for use by the Result_command class."""

416

Source Code for Module lib.frame_order.matrix_ops