maths_fns.rdc

1 ############################################################################### 2 # # 3 # Copyright (C) 2008-2012 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, absolute=False):

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 @keyword absolute: The absolute value or signless RDC flag. 140 @type absolute: int 141 @return: The average RDC value. 142 @rtype: float 143 """ 144 145 # Initial back-calculated RDC value. 146 val = 0.0 147 148 # No weights given. 149 if weights == None: 150 pc = 1.0 / N 151 weights = [pc] * N 152 153 # Missing last weight. 154 if len(weights) < N: 155 pN = 1.0 - sum(weights, axis=0) 156 weights = weights.tolist() 157 weights.append(pN) 158 159 # Back-calculate the RDC. 160 for c in range(N): 161 val = val + weights[c] * dot(vect[c], dot(A, vect[c])) 162 163 # Return the average RDC. 164 if absolute: 165 return abs(dj * val) 166 else: 167 return dj * val

168 169

170 -def ave_rdc_tensor_dDij_dAmn(dj, vect, N, dAi_dAmn, weights=None, absolute=False):

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

235 236

237 -def rdc_tensor(dj, mu, A, absolute=False):

238 """Calculate the RDC, using the 3D alignment tensor. 239 240 The RDC value is:: 241 242 T 243 Dij(theta) = dj . mu_j . Ai . mu_j, 244 245 where: 246 - i is the alignment tensor index, 247 - j is the index over spins, 248 - theta is the parameter vector, 249 - dj is the dipolar constant for spin j, 250 - mu_j i the unit vector corresponding to spin j, 251 - Ai is the alignment tensor. 252 253 The dipolar constant is defined as:: 254 255 dj = 3 / (2pi) d', 256 257 where the factor of 2pi is to convert from units of rad.s^-1 to Hertz, the factor of 3 is 258 associated with the alignment tensor and the pure dipolar constant in SI units is:: 259 260 mu0 gI.gS.h_bar 261 d' = - --- ----------- , 262 4pi r**3 263 264 where: 265 - mu0 is the permeability of free space, 266 - gI and gS are the gyromagnetic ratios of the I and S spins, 267 - h_bar is Dirac's constant which is equal to Planck's constant divided by 2pi, 268 - r is the distance between the two spins. 269 270 271 @param dj: The dipolar constant for spin j. 272 @type dj: float 273 @param mu: The unit XH bond vector. 274 @type mu: numpy rank-1 3D array 275 @param A: The alignment tensor. 276 @type A: numpy rank-2 3D tensor 277 @keyword absolute: The absolute value or signless RDC flag. 278 @type absolute: int 279 @return: The RDC value. 280 @rtype: float 281 """ 282 283 # Return the RDC. 284 if absolute: 285 return abs(dj * dot(mu, dot(A, mu))) 286 else: 287 return dj * dot(mu, dot(A, mu))

288

Source Code for Module maths_fns.rdc