lib.alignment.rdc

1 ############################################################################### 2 # # 3 # Copyright (C) 2008-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 calculation of RDCs.""" 24 25 # Python imports. 26 from numpy import dot, sum 27 28

29 -def ave_rdc_5D(dj, vect, N, A, weights=None):

30 """Calculate the ensemble average RDC, using the 5D tensor. 31 32 This function calculates the average RDC for a set of XH bond vectors from a structural ensemble, using the 5D vector form of the alignment tensor. The formula for this ensemble average RDC value is:: 33 34 _N_ 35 \ 36 Dij(theta) = > pc . RDC_ijc (theta), 37 /__ 38 c=1 39 40 where: 41 - i is the alignment tensor index, 42 - j is the index over spins, 43 - theta is the parameter vector consisting of the alignment tensor parameters {Axx, Ayy, Axy, Axz, Ayz} for each alignment, 44 - c is the index over the states or multiple structures, 45 - N is the total number of states or structures, 46 - pc is the population probability or weight associated with state c (equally weighted to 47 - RDC_ijc is the back-calculated RDC value for alignment tensor i, spin system j and structure c. 48 49 The back-calculated RDC is given by the formula:: 50 51 RDC_ijc(theta) = (x_jc**2 - z_jc**2)Axx_i + (y_jc**2 - z_jc**2)Ayy_i + 2x_jc.y_jc.Axy_i + 2x_jc.z_jc.Axz_i + 2y_jc.z_jc.Ayz_i. 52 53 54 @param dj: The dipolar constant for spin j. 55 @type dj: float 56 @param vect: The unit XH bond vector matrix. The first dimension corresponds to the structural index, the second dimension is the coordinates of the unit vector. 57 @type vect: numpy matrix 58 @param N: The total number of structures. 59 @type N: int 60 @param A: The 5D vector object. The vector format is {Axx, Ayy, Axy, Axz, Ayz}. 61 @type A: numpy 5D vector 62 @param weights: The weights for each member of the ensemble (the last member need not be supplied). 63 @type weights: numpy rank-1 array 64 @return: The average RDC value. 65 @rtype: float 66 """ 67 68 # Initial back-calculated RDC value. 69 val = 0.0 70 71 # No weights given. 72 if weights == None: 73 pc = 1.0 / N 74 weights = [pc] * N 75 76 # Missing last weight. 77 if len(weights) < N: 78 pN = 1.0 - sum(weights, axis=0) 79 weights = weights.tolist() 80 weights.append(pN) 81 82 # Back-calculate the RDC. 83 for c in range(N): 84 val = val + weights[c] * (vect[c, 0]**2 - vect[c, 2]**2)*A[0] + (vect[c, 1]**2 - vect[c, 2]**2)*A[1] + 2.0*vect[c, 0]*vect[c, 1]*A[2] + 2.0*vect[c, 0]*vect[c, 2]*A[3] + 2.0*vect[c, 1]*vect[c, 2]*A[4] 85 86 # Return the average RDC. 87 return dj * val

88 89

90 -def ave_rdc_tensor(dj, vect, N, A, weights=None):

91 """Calculate the ensemble average RDC, using the 3D tensor. 92 93 This function calculates the average RDC for a set of XH bond vectors from a structural ensemble, using the 3D tensorial form of the alignment tensor. The formula for this ensemble average RDC value is:: 94 95 _N_ 96 \ T 97 Dij(theta) = dj > pc . mu_jc . Ai . mu_jc, 98 /__ 99 c=1 100 101 where: 102 - i is the alignment tensor index, 103 - j is the index over spins, 104 - c is the index over the states or multiple structures, 105 - theta is the parameter vector, 106 - dj is the dipolar constant for spin j, 107 - N is the total number of states or structures, 108 - pc is the population probability or weight associated with state c (equally weighted to 1/N if weights are not provided), 109 - mu_jc is the unit vector corresponding to spin j and state c, 110 - Ai is the alignment tensor. 111 112 The dipolar constant is defined as:: 113 114 dj = 3 / (2pi) d', 115 116 where the factor of 2pi is to convert from units of rad.s^-1 to Hertz, the factor of 3 is associated with the alignment tensor and the pure dipolar constant in SI units is:: 117 118 mu0 gI.gS.h_bar 119 d' = - --- ----------- , 120 4pi r**3 121 122 where: 123 - mu0 is the permeability of free space, 124 - gI and gS are the gyromagnetic ratios of the I and S spins, 125 - h_bar is Dirac's constant which is equal to Planck's constant divided by 2pi, 126 - r is the distance between the two spins. 127 128 129 @param dj: The dipolar constant for spin j. 130 @type dj: float 131 @param vect: The unit XH bond vector matrix. The first dimension corresponds to the structural index, the second dimension is the coordinates of the unit vector. 132 @type vect: numpy matrix 133 @param N: The total number of structures. 134 @type N: int 135 @param A: The alignment tensor. 136 @type A: numpy rank-2 3D tensor 137 @keyword weights: The weights for each member of the ensemble (the last member need not be supplied). 138 @type weights: numpy rank-1 array or None 139 @return: The average RDC value. 140 @rtype: float 141 """ 142 143 # Initial back-calculated RDC value. 144 val = 0.0 145 146 # No weights given. 147 if weights == None: 148 pc = 1.0 / N 149 weights = [pc] * N 150 151 # Missing last weight. 152 if len(weights) < N: 153 pN = 1.0 - sum(weights, axis=0) 154 weights = weights.tolist() 155 weights.append(pN) 156 157 # Back-calculate the RDC. 158 for c in range(N): 159 val = val + weights[c] * dot(vect[c], dot(A, vect[c])) 160 161 # Return the average RDC. 162 return dj * val

163 164

165 -def ave_rdc_tensor_dDij_dAmn(dj, vect, N, dAi_dAmn, weights=None):

166 """Calculate the ensemble average RDC gradient element for Amn, using the 3D tensor. 167 168 This function calculates the average RDC gradient for a set of XH bond vectors from a structural ensemble, using the 3D tensorial form of the alignment tensor. The formula for this ensemble average RDC gradient element is:: 169 170 _N_ 171 dDij(theta) \ T dAi 172 ----------- = dj > pc . mu_jc . ---- . mu_jc, 173 dAmn /__ dAmn 174 c=1 175 176 where: 177 - i is the alignment tensor index, 178 - j is the index over spins, 179 - m, the index over the first dimension of the alignment tensor m = {x, y, z}. 180 - n, the index over the second dimension of the alignment tensor n = {x, y, z}, 181 - c is the index over the states or multiple structures, 182 - theta is the parameter vector, 183 - Amn is the matrix element of the alignment tensor, 184 - dj is the dipolar constant for spin j, 185 - N is the total number of states or structures, 186 - pc is the population probability or weight associated with state c (equally weighted to 1/N if weights are not provided), 187 - mu_jc is the unit vector corresponding to spin j and state c, 188 - dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn. 189 190 191 @param dj: The dipolar constant for spin j. 192 @type dj: float 193 @param vect: The unit XH bond vector matrix. The first dimension corresponds to the structural index, the second dimension is the coordinates of the unit vector. 194 @type vect: numpy matrix 195 @param N: The total number of structures. 196 @type N: int 197 @param dAi_dAmn: The alignment tensor derivative with respect to parameter Amn. 198 @type dAi_dAmn: numpy rank-2 3D tensor 199 @keyword weights: The weights for each member of the ensemble (the last member need not be supplied). 200 @type weights: numpy rank-1 array 201 @return: The average RDC gradient element. 202 @rtype: float 203 """ 204 205 # Initial back-calculated RDC gradient. 206 grad = 0.0 207 208 # No weights given. 209 if weights == None: 210 pc = 1.0 / N 211 weights = [pc] * N 212 213 # Missing last weight. 214 if len(weights) < N: 215 pN = 1.0 - sum(weights, axis=0) 216 weights = weights.tolist() 217 weights.append(pN) 218 219 # Back-calculate the RDC gradient element. 220 for c in range(N): 221 grad = grad + weights[c] * dot(vect[c], dot(dAi_dAmn, vect[c])) 222 223 # Return the average RDC gradient element. 224 return dj * grad

225 226

227 -def ave_rdc_tensor_pseudoatom(dj, vect, N, A, weights=None):

228 """Calculate the ensemble and pseudo-atom averaged RDC, using the 3D tensor. 229 230 This function calculates the average RDC for a set of XH bond vectors from a structural ensemble, using the 3D tensorial form of the alignment tensor. The RDC for each pseudo-atom is calculated and then averaged. The formula for this ensemble and pseudo-atom average RDC value is:: 231 232 _N_ _M_ 233 \ 1 \ T 234 Dij(theta) = dj > pc . - > mu_jcd . Ai . mu_jcd, 235 /__ M /__ 236 c=1 d=1 237 238 where: 239 - i is the alignment tensor index, 240 - j is the index over spins, 241 - c is the index over the states or multiple structures, 242 - d is the index over the pseudo-atoms, 243 - theta is the parameter vector, 244 - dj is the dipolar constant for spin j, 245 - N is the total number of states or structures, 246 - M is the total number of pseudo-atoms, 247 - pc is the population probability or weight associated with state c (equally weighted to 1/N if weights are not provided), 248 - mu_jcd is the unit vector corresponding to spin j, state c, and pseudo-atom d, 249 - Ai is the alignment tensor. 250 251 The dipolar constant is defined as:: 252 253 dj = 3 / (2pi) d', 254 255 where the factor of 2pi is to convert from units of rad.s^-1 to Hertz, the factor of 3 is associated with the alignment tensor and the pure dipolar constant in SI units is:: 256 257 mu0 gI.gS.h_bar 258 d' = - --- ----------- , 259 4pi r**3 260 261 where: 262 - mu0 is the permeability of free space, 263 - gI and gS are the gyromagnetic ratios of the I and S spins, 264 - h_bar is Dirac's constant which is equal to Planck's constant divided by 2pi, 265 - r is the distance between the two spins. 266 267 268 @param dj: The dipolar constant for spin j. 269 @type dj: list of float 270 @param vect: The unit vector matrix. The first dimension corresponds to the structural index, the second dimension to the pseudo-atom index, and the third dimension is the coordinates of the unit vector. 271 @type vect: numpy matrix 272 @param N: The total number of structures. 273 @type N: int 274 @param A: The alignment tensor. 275 @type A: numpy rank-2 3D tensor 276 @keyword weights: The weights for each member of the ensemble (the last member need not be supplied). 277 @type weights: numpy rank-1 array or None 278 @return: The average RDC value. 279 @rtype: float 280 """ 281 282 # Calculate M. 283 M = len(dj) 284 285 # Initial back-calculated RDC value. 286 val = 0.0 287 288 # Loop over the pseudo-atoms, calculating the average RDC for the structural ensemble. 289 for d in range(M): 290 val += ave_rdc_tensor(dj[d], vect[:, d,:], N, A, weights) 291 292 # Average. 293 val = val / M 294 295 # Return the averaged value. 296 return val

297 298

299 -def ave_rdc_tensor_pseudoatom_dDij_dAmn(dj, vect, N, dAi_dAmn, weights=None):

300 """Calculate the ensemble and pseudo-atom average RDC gradient element for Amn, using the 3D tensor. 301 302 This function calculates the average RDC gradient for a set of XH bond vectors from a structural ensemble, using the 3D tensorial form of the alignment tensor. The formula for this ensemble average RDC gradient element is:: 303 304 _N_ _M_ 305 dDij(theta) \ 1 \ T dAi 306 ----------- = dj > pc . - > mu_jcd . ---- . mu_jcd, 307 dAmn /__ M /__ dAmn 308 c=1 d=1 309 310 where: 311 - i is the alignment tensor index, 312 - j is the index over spins, 313 - m, the index over the first dimension of the alignment tensor m = {x, y, z}. 314 - n, the index over the second dimension of the alignment tensor n = {x, y, z}, 315 - c is the index over the states or multiple structures, 316 - d is the index over the pseudo-atoms, 317 - theta is the parameter vector, 318 - Amn is the matrix element of the alignment tensor, 319 - dj is the dipolar constant for spin j, 320 - N is the total number of states or structures, 321 - M is the total number of pseudo-atoms, 322 - pc is the population probability or weight associated with state c (equally weighted to 1/N if weights are not provided), 323 - mu_jc is the unit vector corresponding to spin j and state c, 324 - dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn. 325 326 327 @param dj: The dipolar constant for spin j. 328 @type dj: float 329 @param vect: The unit XH bond vector matrix. The first dimension corresponds to the structural index, the second dimension is the coordinates of the unit vector. 330 @type vect: numpy matrix 331 @param N: The total number of structures. 332 @type N: int 333 @param dAi_dAmn: The alignment tensor derivative with respect to parameter Amn. 334 @type dAi_dAmn: numpy rank-2 3D tensor 335 @keyword weights: The weights for each member of the ensemble (the last member need not be supplied). 336 @type weights: numpy rank-1 array 337 @return: The average RDC gradient element. 338 @rtype: float 339 """ 340 341 # Calculate M. 342 M = len(dj) 343 344 # Initial back-calculated RDC value. 345 grad = 0.0 346 347 # Loop over the pseudo-atoms, calculating the average RDC for the structural ensemble. 348 for d in range(M): 349 grad += ave_rdc_tensor_dDij_dAmn(dj[d], vect[:, d,:], N, dAi_dAmn, weights) 350 351 # Average. 352 grad = grad / M 353 354 # Return the average RDC gradient element. 355 return grad

356 357

358 -def rdc_tensor(dj, mu, A):

359 """Calculate the RDC, using the 3D alignment tensor. 360 361 The RDC value is:: 362 363 T 364 Dij(theta) = dj . mu_j . Ai . mu_j, 365 366 where: 367 - i is the alignment tensor index, 368 - j is the index over spins, 369 - theta is the parameter vector, 370 - dj is the dipolar constant for spin j, 371 - mu_j i the unit vector corresponding to spin j, 372 - Ai is the alignment tensor. 373 374 The dipolar constant is defined as:: 375 376 dj = 3 / (2pi) d', 377 378 where the factor of 2pi is to convert from units of rad.s^-1 to Hertz, the factor of 3 is 379 associated with the alignment tensor and the pure dipolar constant in SI units is:: 380 381 mu0 gI.gS.h_bar 382 d' = - --- ----------- , 383 4pi r**3 384 385 where: 386 - mu0 is the permeability of free space, 387 - gI and gS are the gyromagnetic ratios of the I and S spins, 388 - h_bar is Dirac's constant which is equal to Planck's constant divided by 2pi, 389 - r is the distance between the two spins. 390 391 392 @param dj: The dipolar constant for spin j. 393 @type dj: float 394 @param mu: The unit XH bond vector. 395 @type mu: numpy rank-1 3D array 396 @param A: The alignment tensor. 397 @type A: numpy rank-2 3D tensor 398 @return: The RDC value. 399 @rtype: float 400 """ 401 402 # Return the RDC. 403 return dj * dot(mu, dot(A, mu))

404

Source Code for Module lib.alignment.rdc