maths_fns.rdc

1 ############################################################################### 2 # # 3 # Copyright (C) 2008, 2010 Edward d'Auvergne # 4 # # 5 # This file is part of the program relax. # 6 # # 7 # relax 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 2 of the License, or # 10 # (at your option) any later version. # 11 # # 12 # relax 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 relax; if not, write to the Free Software # 19 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA # 20 # # 21 ############################################################################### 22 23 # Module docstring. 24 """Module for the calculation of RDCs.""" 25 26 # Python imports. 27 from numpy import dot, sum 28 29

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

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

89 90

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

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

164 165

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

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

226 227

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

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

274

Source Code for Module maths_fns.rdc