lib.alignment.pcs

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 pseudocontact shifts.""" 24 25 # Python imports. 26 from math import pi 27 from numpy import dot, sum 28 29 # relax module imports. 30 from lib.physical_constants import kB, mu0 31 32

33 -def ave_pcs_tensor(dj, vect, N, A, weights=None):

34 """Calculate the ensemble average PCS, using the 3D tensor. 35 36 This function calculates the average PCS 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 PCS value is:: 37 38 _N_ 39 \ T 40 <delta_ij(theta)> = > pc . djc . mu_jc . Ai . mu_jc, 41 /__ 42 c=1 43 44 where: 45 - i is the alignment tensor index, 46 - j is the index over spins, 47 - c is the index over the states or multiple structures, 48 - N is the total number of states or structures, 49 - theta is the parameter vector, 50 - djc is the PCS constant for spin j and state c, 51 - pc is the population probability or weight associated with state c (equally weighted to 1/N if weights are not provided), 52 - mu_jc is the unit vector corresponding to spin j and state c, 53 - Ai is the alignment tensor. 54 55 The PCS constant is defined as:: 56 57 mu0 15kT 1 58 dj = --- ----- ---- , 59 4pi Bo**2 r**3 60 61 where: 62 - mu0 is the permeability of free space, 63 - k is Boltzmann's constant, 64 - T is the absolute temperature, 65 - Bo is the magnetic field strength, 66 - r is the distance between the paramagnetic centre (electron spin) and the nuclear spin. 67 68 69 @param dj: The PCS constants for each structure c for spin j. This should be an array with indices corresponding to c. 70 @type dj: numpy rank-1 array 71 @param vect: The electron-nuclear unit vector matrix. The first dimension corresponds to the structural index, the second dimension is the coordinates of the unit vector. The vectors should be parallel to the vector connecting the paramagnetic centre to the nuclear spin. 72 @type vect: numpy matrix 73 @param N: The total number of structures. 74 @type N: int 75 @param A: The alignment tensor. 76 @type A: numpy rank-2 3D tensor 77 @param weights: The weights for each member of the ensemble (the last member need not be supplied). 78 @type weights: numpy rank-1 array 79 @return: The average PCS value. 80 @rtype: float 81 """ 82 83 # Initial back-calculated PCS value. 84 val = 0.0 85 86 # No weights given. 87 if weights == None: 88 pc = 1.0 / N 89 weights = [pc] * N 90 91 # Missing last weight. 92 if len(weights) < N: 93 pN = 1.0 - sum(weights, axis=0) 94 weights = weights.tolist() 95 weights.append(pN) 96 97 # Back-calculate the PCS. 98 for c in range(N): 99 val = val + weights[c] * dj[c] * dot(vect[c], dot(A, vect[c])) 100 101 # Return the average PCS. 102 return val

103 104

105 -def ave_pcs_tensor_ddeltaij_dAmn(dj, vect, N, dAi_dAmn, weights=None):

106 """Calculate the ensemble average PCS gradient element for Amn, using the 3D tensor. 107 108 This function calculates the alignment tensor parameter part of the average PCS gradient for a set of electron-nuclear spin unit vectors (paramagnetic to the nuclear spin) from a structural ensemble, using the 3D tensorial form of the alignment tensor. The formula for this ensemble average PCS gradient element is:: 109 110 _N_ 111 ddelta_ij(theta) \ T dAi 112 ---------------- = > pc . djc . mu_jc . ---- . mu_jc, 113 dAmn /__ dAmn 114 c=1 115 116 where: 117 - i is the alignment tensor index, 118 - j is the index over spins, 119 - m, the index over the first dimension of the alignment tensor m = {x, y, z}. 120 - n, the index over the second dimension of the alignment tensor n = {x, y, z}, 121 - c is the index over the states or multiple structures, 122 - theta is the parameter vector, 123 - Amn is the matrix element of the alignment tensor, 124 - djc is the PCS constant for spin j and state c, 125 - N is the total number of states or structures, 126 - pc is the population probability or weight associated with state c (equally weighted to 1/N if weights are not provided), 127 - mu_jc is the unit vector corresponding to spin j and state c, 128 - dAi/dAmn is the partial derivative of the alignment tensor with respect to element Amn. 129 130 131 @param dj: The PCS constants for each structure c for spin j. This should be an array with indices corresponding to c. 132 @type dj: numpy rank-1 array 133 @param vect: The electron-nuclear unit vector matrix. The first dimension corresponds to the structural index, the second dimension is the coordinates of the unit vector. The vectors should be parallel to the vector connecting the paramagnetic centre to the nuclear spin. 134 @type vect: numpy matrix 135 @param N: The total number of structures. 136 @type N: int 137 @param dAi_dAmn: The alignment tensor derivative with respect to parameter Amn. 138 @type dAi_dAmn: numpy rank-2 3D tensor 139 @param weights: The weights for each member of the ensemble (the last member need not be supplied). 140 @type weights: numpy rank-1 array 141 @return: The average PCS gradient element. 142 @rtype: float 143 """ 144 145 # Initial back-calculated PCS gradient. 146 grad = 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 PCS gradient element. 160 for c in range(N): 161 grad = grad + weights[c] * dj[c] * dot(vect[c], dot(dAi_dAmn, vect[c])) 162 163 # Return the average PCS gradient element. 164 return grad

165 166

167 -def ave_pcs_tensor_ddeltaij_dc(ddj=None, dj=None, r=None, unit_vect=None, N=None, Ai=None, dr_dc=None, weights=None):

168 """Calculate the ensemble average PCS gradient element for the paramagnetic centre coordinate c, using the 3D tensor. 169 170 This function calculates the paramagnetic centre coordinate part of the average PCS gradient for a set of electron-nuclear spin unit vectors (paramagnetic to the nuclear spin) from a structural ensemble, using the 3D tensorial form of the alignment tensor. The formula for this ensemble average PCS gradient element is:: 171 172 _N_ 173 ddelta_ij(theta) \ / ddjc dr_jcT dr_jc \ 174 ---------------- = > pc . | ----.r_jcT.Ai.r_jc + djc.------.Ai.r_jc + djc.r_jcT.Ai.----- | , 175 dxi /__ \ dxi dxi dxi / 176 c=1 177 178 where the last two terms in the sum are equal due to the symmetry of the alignment tensor, and: 179 - xi are the paramagnetic position coordinates {x0, x1, x2}, 180 - i is the alignment tensor index, 181 - j is the index over spins, 182 - c is the index over the states or multiple structures, 183 - theta is the parameter vector, 184 - djc is the PCS constant for spin j and state c, 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 - r_jc is the vector corresponding to spin j and state c, 188 189 and where:: 190 191 ddjc mu0 15kT 5 (si - xi) 192 ---- = --- ----- --------------------------------------------- , 193 dxi 4pi Bo**2 ((sx-x0)**2 + (sy-x1)**2 + (sz-x2)**2)**(7/2) 194 195 where {sx, sy, sz} are the spin atomic coordinates, and:: 196 197 dr | 1 | dr | 0 | dr | 0 | 198 --- = - | 0 | , --- = - | 1 | , --- = - | 0 | . 199 dx0 | 0 | dx1 | 0 | dx2 | 1 | 200 201 The pseudocontact shift constant is defined here as:: 202 203 mu0 15kT 1 204 djc = --- ----- ------ , 205 4pi Bo**2 rjc**5 206 207 208 @keyword ddj: The PCS constant gradient for each structure c for spin j. This should be an array with indices corresponding to c. 209 @type ddj: numpy rank-1 array 210 @keyword dj: The PCS constants for each structure c for spin j. This should be an array with indices corresponding to c. 211 @type dj: numpy rank-1 array 212 @keyword r: The distance between the paramagnetic centre and the spin (in meters). 213 @type r: float 214 @keyword vect: The electron-nuclear unit vector matrix. The first dimension corresponds to the structural index, the second dimension is the coordinates of the unit vector. The vectors should be parallel to the vector connecting the paramagnetic centre to the nuclear spin. 215 @type vect: numpy matrix 216 @keyword N: The total number of structures. 217 @type N: int 218 @keyword Ai: The alignment tensor i. 219 @type Ai: numpy rank-2, 3D tensor 220 @keyword weights: The weights for each member of the ensemble (the last member need not be supplied). 221 @type weights: numpy rank-1 array 222 @return: The average PCS gradient element. 223 @rtype: float 224 """ 225 226 # Initial back-calculated PCS gradient. 227 grad = 0.0 228 229 # No weights given. 230 if weights == None: 231 pc = 1.0 / N 232 weights = [pc] * N 233 234 # Missing last weight. 235 if len(weights) < N: 236 pN = 1.0 - sum(weights, axis=0) 237 weights = weights.tolist() 238 weights.append(pN) 239 240 # Loop over each state. 241 for c in range(N): 242 # Recreate the full length vector. 243 vect = unit_vect[c] * r[c] 244 245 # Modified PCS constant (from r**-3 to r**-5). 246 const = dj[c] / r[c]**2 247 248 # Back-calculate the PCS gradient element. 249 grad += weights[c] * (ddj[c] * dot(vect, dot(Ai, vect)) + 2.0 * const * dot(dr_dc, dot(Ai, vect))) 250 251 # Return the average PCS gradient element, converted back to the Angstrom scale at the coordinates are in Angstrom units. 252 return grad * 1e-10

253 254

255 -def pcs_constant_grad(T=None, Bo=None, r=None, unit_vect=None, grad=None):

256 """Calculate the PCS constant gradient with respect to the paramagnetic centre position. 257 258 The pseudocontact shift constant is defined here as:: 259 260 mu0 15kT 1 261 djc = --- ----- ------ , 262 4pi Bo**2 rjc**5 263 264 where: 265 - mu0 is the permeability of free space, 266 - k is Boltzmann's constant, 267 - T is the absolute temperature, 268 - Bo is the magnetic field strength, 269 - r is the distance between the paramagnetic centre (electron spin) and the nuclear spin. 270 271 The 5th power of the distance is used to simplify the PCS derivative. The pseudocontact shift constant derivative is:: 272 273 ddjc mu0 15kT 5 (si - xi) 274 ---- = --- ----- --------------------------------------------- , 275 dxi 4pi Bo**2 ((sx-x0)**2 + (sy-x1)**2 + (sz-x2)**2)**(7/2) 276 277 where: 278 - {x0, x1, x2} are the paramagnetic centre coordinates, 279 - {sx, sy, sz} are the spin atomic coordinates. 280 281 282 @keyword T: The temperature in kelvin. 283 @type T: float 284 @keyword Bo: The magnetic field strength. 285 @type Bo: float 286 @keyword r: The distance between the paramagnetic centre and the spin (in meters). 287 @type r: float 288 @keyword unit_vect: The paramagnetic centre to spin unit vector. 289 @type unit_vect: numpy rank-1, 3D array 290 @keyword grad: The gradient component to update. The indices {0, 1, 2} match the {dx, dy, dz} derivatives. 291 @type grad: numpy rank-1, 3D array 292 """ 293 294 # Recreate the full length vector. 295 vect = unit_vect * r 296 297 # Calculate the invariant part. 298 a = 18.75 * mu0 / pi * kB * T / Bo**2 299 300 # Calculate the coordinate part. 301 b = r**7 302 303 # Combine. 304 for i in range(3): 305 grad[i] = a * vect[i] / b

306 307

308 -def pcs_tensor(dj, mu, A):

309 """Calculate the PCS, using the 3D alignment tensor. 310 311 The PCS value is:: 312 313 T 314 delta_ij(theta) = dj . mu_j . Ai . mu_j, 315 316 where: 317 - i is the alignment tensor index, 318 - j is the index over spins, 319 - theta is the parameter vector, 320 - dj is the PCS constant for spin j, 321 - mu_j i the unit vector corresponding to spin j, 322 - Ai is the alignment tensor. 323 324 The PCS constant is defined as:: 325 326 mu0 15kT 1 327 dj = --- ----- ---- , 328 4pi Bo**2 r**3 329 330 where: 331 - mu0 is the permeability of free space, 332 - k is Boltzmann's constant, 333 - T is the absolute temperature, 334 - Bo is the magnetic field strength, 335 - r is the distance between the paramagnetic centre (electron spin) and the nuclear spin. 336 337 338 @param dj: The PCS constant for spin j. 339 @type dj: float 340 @param mu: The unit vector connecting the electron and nuclear spins. 341 @type mu: numpy rank-1 3D array 342 @param A: The alignment tensor. 343 @type A: numpy rank-2 3D tensor 344 @return: The PCS value. 345 @rtype: float 346 """ 347 348 # Return the PCS. 349 return dj * dot(mu, dot(A, mu))

350

Source Code for Module lib.alignment.pcs