lib.frame_order.double

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 double rotor frame order model.""" 24 25 # Python module imports. 26 from math import cos, pi, sin, sqrt 27 from numpy import dot, inner, sinc, transpose 28 29 # relax module imports. 30 from lib.frame_order.matrix_ops import rotate_daeg 31 32

33 -def compile_2nd_matrix_double_rotor(matrix, Rx2_eigen, smax, smaxb):

34 """Generate the rotated 2nd degree Frame Order matrix for the double rotor model. 35 36 The cone axis is assumed to be parallel to the z-axis in the eigenframe. 37 38 39 @param matrix: The Frame Order matrix, 2nd degree to be populated. 40 @type matrix: numpy 9D, rank-2 array 41 @param Rx2_eigen: The Kronecker product of the eigenframe rotation matrix with itself. 42 @type Rx2_eigen: numpy 9D, rank-2 array 43 @param smax: The maximum torsion angle for the first rotor. 44 @type smax: float 45 @param smaxb: The maximum torsion angle for the second rotor. 46 @type smaxb: float 47 """ 48 49 # Zeros. 50 for i in range(9): 51 for j in range(9): 52 matrix[i, j] = 0.0 53 54 # Repetitive trig calculations. 55 sinc_smax = sinc(smax/pi) 56 sinc_2smax = sinc(2.0*smax/pi) 57 58 # Diagonal. 59 matrix[0, 0] = (sinc_2smax + 1.0) / 2.0 60 matrix[1, 1] = matrix[0, 0] 61 matrix[2, 2] = sinc_smax 62 matrix[3, 3] = matrix[0, 0] 63 matrix[4, 4] = matrix[0, 0] 64 matrix[5, 5] = matrix[2, 2] 65 matrix[6, 6] = matrix[2, 2] 66 matrix[7, 7] = matrix[2, 2] 67 matrix[8, 8] = 1.0 68 69 # Off diagonal set 1. 70 matrix[0, 4] = matrix[4, 0] = -(sinc_2smax - 1.0) / 2.0 71 72 # Off diagonal set 2. 73 matrix[1, 3] = matrix[3, 1] = -matrix[0, 4] 74 75 # Rotate and return the frame order matrix. 76 return rotate_daeg(matrix, Rx2_eigen)

77 78

79 -def pcs_numeric_int_double_rotor(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, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None, pcs_theta=None, pcs_theta_err=None, missing_pcs=None, error_flag=False):

80 """The averaged PCS value via numerical integration for the double rotor frame order model. 81 82 @keyword points: The Sobol points in the torsion-tilt angle space. 83 @type points: numpy rank-2, 3D array 84 @keyword sigma_max: The maximum opening angle for the first rotor. 85 @type sigma_max: float 86 @keyword sigma_max_2: The maximum opening angle for the second rotor. 87 @type sigma_max_2: float 88 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units). 89 @type c: numpy rank-1 array 90 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 91 @type full_in_ref_frame: numpy rank-1 array 92 @keyword r_pivot_atom: The pivot point to atom vector. 93 @type r_pivot_atom: numpy rank-2, 3D array 94 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector. 95 @type r_pivot_atom_rev: numpy rank-2, 3D array 96 @keyword r_ln_pivot: The lanthanide position to pivot point vector. 97 @type r_ln_pivot: numpy rank-2, 3D array 98 @keyword A: The full alignment tensor of the non-moving domain. 99 @type A: numpy rank-2, 3D array 100 @keyword R_eigen: The eigenframe rotation matrix. 101 @type R_eigen: numpy rank-2, 3D array 102 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations). 103 @type RT_eigen: numpy rank-2, 3D array 104 @keyword Ri_prime: The empty rotation matrix for the in-frame rotor motion, used to calculate the PCS for each state i in the numerical integration. 105 @type Ri_prime: numpy rank-2, 3D array 106 @keyword pcs_theta: The storage structure for the back-calculated PCS values. 107 @type pcs_theta: numpy rank-2 array 108 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors. 109 @type pcs_theta_err: numpy rank-2 array 110 @keyword missing_pcs: A structure used to indicate which PCS values are missing. 111 @type missing_pcs: numpy rank-2 array 112 @keyword error_flag: A flag which if True will cause the PCS errors to be estimated and stored in pcs_theta_err. 113 @type error_flag: bool 114 """ 115 116 # Clear the data structures. 117 for i in range(len(pcs_theta)): 118 for j in range(len(pcs_theta[i])): 119 pcs_theta[i, j] = 0.0 120 pcs_theta_err[i, j] = 0.0 121 122 # Loop over the samples. 123 num = 0 124 for i in range(len(points)): 125 # Unpack the point. 126 sigma, sigma2 = points[i] 127 128 # Outside of the distribution, so skip the point. 129 if sigma > sigma_max or sigma < -sigma_max: 130 continue 131 132 # Calculate the PCSs for this state. 133 pcs_pivot_motion_double_rotor(sigma_i=sigma, 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) 134 135 # Increment the number of points. 136 num += 1 137 138 # Calculate the PCS and error. 139 for i in range(len(pcs_theta)): 140 for j in range(len(pcs_theta[i])): 141 # The average PCS. 142 pcs_theta[i, j] = c[i] * pcs_theta[i, j] / float(num) 143 144 # The error. 145 if error_flag: 146 pcs_theta_err[i, j] = abs(pcs_theta_err[i, j] / float(num) - pcs_theta[i, j]**2) / float(num) 147 pcs_theta_err[i, j] = c[i] * sqrt(pcs_theta_err[i, j]) 148 print("%8.3f +/- %-8.3f" % (pcs_theta[i, j]*1e6, pcs_theta_err[i, j]*1e6))

149 150

151 -def pcs_pivot_motion_double_rotor(sigma_i=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):

152 """Calculate the PCS value after a pivoted motion for the double rotor model. 153 154 @keyword sigma_i: The rotor angle for state i. 155 @type sigma_i: float 156 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 157 @type full_in_ref_frame: numpy rank-1 array 158 @keyword r_pivot_atom: The pivot point to atom vector. 159 @type r_pivot_atom: numpy rank-2, 3D array 160 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector. 161 @type r_pivot_atom_rev: numpy rank-2, 3D array 162 @keyword r_ln_pivot: The lanthanide position to pivot point vector. 163 @type r_ln_pivot: numpy rank-2, 3D array 164 @keyword A: The full alignment tensor of the non-moving domain. 165 @type A: numpy rank-2, 3D array 166 @keyword R_eigen: The eigenframe rotation matrix. 167 @type R_eigen: numpy rank-2, 3D array 168 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations). 169 @type RT_eigen: numpy rank-2, 3D array 170 @keyword Ri_prime: The empty rotation matrix for the in-frame rotor motion for state i. 171 @type Ri_prime: numpy rank-2, 3D array 172 @keyword pcs_theta: The storage structure for the back-calculated PCS values. 173 @type pcs_theta: numpy rank-2 array 174 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors. 175 @type pcs_theta_err: numpy rank-2 array 176 @keyword missing_pcs: A structure used to indicate which PCS values are missing. 177 @type missing_pcs: numpy rank-2 array 178 @keyword error_flag: A flag which if True will cause the PCS errors to be estimated and stored in pcs_theta_err. 179 @type error_flag: bool 180 """ 181 182 # The rotation matrix. 183 c_sigma = cos(sigma_i) 184 s_sigma = sin(sigma_i) 185 Ri_prime[0, 0] = c_sigma 186 Ri_prime[0, 1] = -s_sigma 187 Ri_prime[0, 2] = 0.0 188 Ri_prime[1, 0] = s_sigma 189 Ri_prime[1, 1] = c_sigma 190 Ri_prime[1, 2] = 0.0 191 Ri_prime[2, 0] = 0.0 192 Ri_prime[2, 1] = 0.0 193 Ri_prime[2, 2] = 1.0 194 195 # The rotation. 196 R_i = dot(R_eigen, dot(Ri_prime, RT_eigen)) 197 198 # Pre-calculate all the new vectors (forwards and reverse). 199 rot_vect_rev = transpose(dot(R_i, r_pivot_atom_rev) + r_ln_pivot) 200 rot_vect = transpose(dot(R_i, r_pivot_atom) + r_ln_pivot) 201 202 # Loop over the atoms. 203 for j in range(len(r_pivot_atom[0])): 204 # The vector length (to the 5th power). 205 length_rev = 1.0 / sqrt(inner(rot_vect_rev[j], rot_vect_rev[j]))**5 206 length = 1.0 / sqrt(inner(rot_vect[j], rot_vect[j]))**5 207 208 # Loop over the alignments. 209 for i in range(len(pcs_theta)): 210 # Skip missing data. 211 if missing_pcs[i, j]: 212 continue 213 214 # The projection. 215 if full_in_ref_frame[i]: 216 proj = dot(rot_vect[j], dot(A[i], rot_vect[j])) 217 length_i = length 218 else: 219 proj = dot(rot_vect_rev[j], dot(A[i], rot_vect_rev[j])) 220 length_i = length_rev 221 222 # The PCS. 223 pcs_theta[i, j] += proj * length_i 224 225 # The square. 226 if error_flag: 227 pcs_theta_err[i, j] += (proj * length_i)**2

228

Source Code for Module lib.frame_order.double_rotor