lib.dispersion.ns_mmq

1 ############################################################################### 2 # # 3 # Copyright (C) 2013 Mathilde Lescanne # 4 # Copyright (C) 2013 Dominique Marion # 5 # Copyright (C) 2013-2014,2019 Edward d'Auvergne # 6 # Copyright (C) 2014 Troels E. Linnet # 7 # # 8 # This file is part of the program relax (http://www.nmr-relax.com). # 9 # # 10 # This program is free software: you can redistribute it and/or modify # 11 # it under the terms of the GNU General Public License as published by # 12 # the Free Software Foundation, either version 3 of the License, or # 13 # (at your option) any later version. # 14 # # 15 # This program is distributed in the hope that it will be useful, # 16 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 17 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 18 # GNU General Public License for more details. # 19 # # 20 # You should have received a copy of the GNU General Public License # 21 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 22 # # 23 ############################################################################### 24 25 # Module docstring. 26 """The numeric solution for the 3-site Bloch-McConnell equations for MMQ CPMG data, the U{NS MMQ 3-site linear<http://wiki.nmr-relax.com/NS_MMQ_3-site_linear>} and U{NS MMQ 3-site<http://wiki.nmr-relax.com/NS_MMQ_3-site>} models. 27 28 Description 29 =========== 30 31 This handles proton-heteronuclear SQ, ZQ, DQ and MQ CPMG data. 32 33 34 References 35 ========== 36 37 It uses the equations of: 38 39 - Dmitry M. Korzhnev, Philipp Neudecker, Anthony Mittermaier, Vladislav Yu. Orekhov, and Lewis E. Kay (2005). Multiple-site exchange in proteins studied with a suite of six NMR relaxation dispersion experiments: An application to the folding of a Fyn SH3 domain mutant. I{J. Am. Chem. Soc.}, B{127}, 15602-15611. (U{DOI: 10.1021/ja054550e<http://dx.doi.org/10.1021/ja054550e>}). 40 41 42 Links 43 ===== 44 45 More information on the NS MMQ 3-site linear model can be found in the: 46 47 - U{relax wiki<http://wiki.nmr-relax.com/NS_MMQ_3-site_linear>}, 48 - U{relax manual<http://www.nmr-relax.com/manual/The_NS_MMQ_3_site_linear_model.html>}, 49 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_MMQ_3-site_linear>}. 50 51 More information on the NS MMQ 3-site model can be found in the: 52 53 - U{relax wiki<http://wiki.nmr-relax.com/NS_MMQ_3-site>}, 54 - U{relax manual<http://www.nmr-relax.com/manual/The_NS_MMQ_3_site_model.html>}, 55 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_MMQ_3-site>}. 56 """ 57 58 # Python module imports. 59 from math import floor 60 from numpy import array, conj, dot, einsum, float64, log, multiply 61 from numpy.linalg import matrix_power 62 63 # relax module imports. 64 from lib.float import isNaN 65 from lib.dispersion.matrix_exponential import matrix_exponential 66 67 # Repetitive calculations (to speed up calculations). 68 # R20. 69 m_r20a = array([ 70 [-1, 0, 0], 71 [ 0, 0, 0], 72 [ 0, 0, 0]], float64) 73 74 m_r20b = array([ 75 [ 0, 0, 0], 76 [ 0, -1, 0], 77 [ 0, 0, 0]], float64) 78 79 m_r20c = array([ 80 [ 0, 0, 0], 81 [ 0, 0, 0], 82 [ 0, 0, -1]], float64) 83 84 # dw. 85 m_dw_AB = array([ 86 [ 0, 0, 0], 87 [ 0, 1, 0], 88 [ 0, 0, 0]], float64) 89 90 m_dw_AC = array([ 91 [ 0, 0, 0], 92 [ 0, 0, 0], 93 [ 0, 0, 1]], float64) 94 95 # k_x. 96 m_k_AB = array([ 97 [-1, 0, 0], 98 [ 1, 0, 0], 99 [ 0, 0, 0]], float64) 100 101 m_k_BA = array([ 102 [ 0, 1, 0], 103 [ 0, -1, 0], 104 [ 0, 0, 0]], float64) 105 106 m_k_BC = array([ 107 [ 0, 0, 0], 108 [ 0, -1, 0], 109 [ 0, 1, 0]], float64) 110 111 m_k_CB = array([ 112 [ 0, 0, 0], 113 [ 0, 0, 1], 114 [ 0, 0, -1]], float64) 115 116 m_k_AC = array([ 117 [-1, 0, 0], 118 [ 0, 0, 0], 119 [ 1, 0, 0]], float64) 120 121 m_k_CA = array([ 122 [ 0, 0, 1], 123 [ 0, 0, 0], 124 [ 0, 0, -1]], float64) 125 126

127 -def rmmq_3site_rankN(R20A=None, R20B=None, R20C=None, dw_AB=None, dw_AC=None, k_AB=None, k_BA=None, k_BC=None, k_CB=None, k_AC=None, k_CA=None, tcp=None):

128 """The Bloch-McConnell matrix for 3-site exchange. 129 130 @keyword R20A: The transverse, spin-spin relaxation rate for state A. 131 @type R20A: numpy float array of rank [NS][NM][NO][ND] 132 @keyword R20B: The transverse, spin-spin relaxation rate for state B. 133 @type R20B: numpy float array of rank [NS][NM][NO][ND] 134 @keyword R20C: The transverse, spin-spin relaxation rate for state C. 135 @type R20C: numpy float array of rank [NS][NM][NO][ND] 136 @keyword dw_AB: The combined chemical exchange difference parameters between states A and B in rad/s. This can be any combination of dw and dwH. 137 @type dw_AB: numpy float array of rank [NS][NM][NO][ND] 138 @keyword dw_AC: The combined chemical exchange difference parameters between states A and C in rad/s. This can be any combination of dw and dwH. 139 @type dw_AC: numpy float array of rank [NS][NM][NO][ND] 140 @keyword k_AB: The rate of exchange from site A to B (rad/s). 141 @type k_AB: float 142 @keyword k_BA: The rate of exchange from site B to A (rad/s). 143 @type k_BA: float 144 @keyword k_BC: The rate of exchange from site B to C (rad/s). 145 @type k_BC: float 146 @keyword k_CB: The rate of exchange from site C to B (rad/s). 147 @type k_CB: float 148 @keyword k_AC: The rate of exchange from site A to C (rad/s). 149 @type k_AC: float 150 @keyword k_CA: The rate of exchange from site C to A (rad/s). 151 @type k_CA: float 152 @keyword tcp: The tau_CPMG times (1 / 4.nu1). 153 @type tcp: numpy float array of rank [NE][NS][NM][NO][ND] 154 """ 155 156 # The first row. 157 #matrix[0, 0] = -k_AB - k_AC - R20A 158 #matrix[0, 1] = k_BA 159 #matrix[0, 2] = k_CA 160 161 # The second row. 162 #matrix[1, 0] = k_AB 163 #matrix[1, 1] = -k_BA - k_BC + 1.j*dw_AB - R20B 164 #matrix[1, 2] = k_CB 165 166 # The third row. 167 #matrix[2, 0] = k_AC 168 #matrix[2, 1] = k_BC 169 #matrix[2, 2] = -k_CB - k_CA + 1.j*dw_AC - R20C 170 171 # Pre-multiply with tcp. 172 r20a_tcp = R20A * tcp 173 r20b_tcp = R20B * tcp 174 r20c_tcp = R20C * tcp 175 176 # Complex dw. 177 dw_AB_C_tcp = dw_AB * tcp * 1j 178 dw_AC_C_tcp = dw_AC * tcp * 1j 179 180 k_AB_tcp = k_AB * tcp 181 k_BA_tcp = k_BA * tcp 182 k_BC_tcp = k_BC * tcp 183 k_CB_tcp = k_CB * tcp 184 k_AC_tcp = k_AC * tcp 185 k_CA_tcp = k_CA * tcp 186 187 # Multiply and expand. 188 m_r20a_tcp = multiply.outer( r20a_tcp, m_r20a ) 189 m_r20b_tcp = multiply.outer( r20b_tcp, m_r20b ) 190 m_r20c_tcp = multiply.outer( r20c_tcp, m_r20c ) 191 192 # Multiply and expand. 193 m_dw_AB_C_tcp = multiply.outer( dw_AB_C_tcp, m_dw_AB ) 194 m_dw_AC_C_tcp = multiply.outer( dw_AC_C_tcp, m_dw_AC ) 195 196 # Multiply and expand. 197 m_k_AB_tcp = multiply.outer( k_AB_tcp, m_k_AB ) 198 m_k_BA_tcp = multiply.outer( k_BA_tcp, m_k_BA ) 199 m_k_BC_tcp = multiply.outer( k_BC_tcp, m_k_BC ) 200 m_k_CB_tcp = multiply.outer( k_CB_tcp, m_k_CB ) 201 m_k_AC_tcp = multiply.outer( k_AC_tcp, m_k_AC ) 202 m_k_CA_tcp = multiply.outer( k_CA_tcp, m_k_CA ) 203 204 # Collect matrix. 205 matrix = (m_r20a_tcp + m_r20b_tcp + m_r20c_tcp 206 + m_dw_AB_C_tcp + m_dw_AC_C_tcp 207 + m_k_AB_tcp + m_k_BA_tcp + m_k_BC_tcp 208 + m_k_CB_tcp + m_k_AC_tcp + m_k_CA_tcp) 209 210 return matrix

211 212

213 -def r2eff_ns_mmq_3site_mq(M0=None, F_vector=array([1, 0, 0], float64), R20A=None, R20B=None, R20C=None, pA=None, pB=None, dw_AB=None, dw_BC=None, dwH_AB=None, dwH_BC=None, kex_AB=None, kex_BC=None, kex_AC=None, inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None):

214 """The 3-site numerical solution to the Bloch-McConnell equation for MQ data. 215 216 The notation used here comes from: 217 218 - Dmitry M. Korzhnev, Philipp Neudecker, Anthony Mittermaier, Vladislav Yu. Orekhov, and Lewis E. Kay (2005). Multiple-site exchange in proteins studied with a suite of six NMR relaxation dispersion experiments: An application to the folding of a Fyn SH3 domain mutant. J. Am. Chem. Soc., 127, 15602-15611. (doi: http://dx.doi.org/10.1021/ja054550e). 219 220 and: 221 222 - Dmitry M. Korzhnev, Philipp Neudecker, Anthony Mittermaier, Vladislav Yu. Orekhov, and Lewis E. Kay (2005). Multiple-site exchange in proteins studied with a suite of six NMR relaxation dispersion experiments: An application to the folding of a Fyn SH3 domain mutant. J. Am. Chem. Soc., 127, 15602-15611. (doi: http://dx.doi.org/10.1021/ja054550e). 223 224 This function calculates and stores the R2eff values. 225 226 227 @keyword M0: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 228 @type M0: numpy float64, rank-1, 7D array 229 @keyword F_vector: The observable magnitisation vector. This defaults to [1, 0] for X observable magnitisation. 230 @type F_vector: numpy rank-1, 3D float64 array 231 @keyword R20A: The transverse, spin-spin relaxation rate for state A. 232 @type R20A: numpy float array of rank [NS][NM][NO][ND] 233 @keyword R20B: The transverse, spin-spin relaxation rate for state B. 234 @type R20B: numpy float array of rank [NS][NM][NO][ND] 235 @keyword R20C: The transverse, spin-spin relaxation rate for state C. 236 @type R20C: numpy float array of rank [NS][NM][NO][ND] 237 @keyword pA: The population of state A. 238 @type pA: float 239 @keyword pB: The population of state B. 240 @type pB: float 241 @keyword dw_AB: The chemical exchange difference between states A and B in rad/s. 242 @type dw_AB: numpy float array of rank [NS][NM][NO][ND] 243 @keyword dw_BC: The chemical exchange difference between states B and C in rad/s. 244 @type dw_BC: numpy float array of rank [NS][NM][NO][ND] 245 @keyword dwH_AB: The proton chemical exchange difference between states A and B in rad/s. 246 @type dwH_AB: numpy float array of rank [NS][NM][NO][ND] 247 @keyword dwH_BC: The proton chemical exchange difference between states B and C in rad/s. 248 @type dwH_BC: numpy float array of rank [NS][NM][NO][ND] 249 @keyword kex_AB: The exchange rate between sites A and B for 3-site exchange with kex_AB = k_AB + k_BA (rad.s^-1) 250 @type kex_AB: float 251 @keyword kex_BC: The exchange rate between sites A and C for 3-site exchange with kex_AC = k_AC + k_CA (rad.s^-1) 252 @type kex_BC: float 253 @keyword kex_AC: The exchange rate between sites B and C for 3-site exchange with kex_BC = k_BC + k_CB (rad.s^-1) 254 @type kex_AC: float 255 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds). 256 @type inv_tcpmg: float 257 @keyword tcp: The tau_CPMG times (1 / 4.nu1). 258 @type tcp: numpy float array of rank [NS][NM][NO][ND] 259 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies. 260 @type back_calc: numpy float array of rank [NS][NM][NO][ND] 261 @keyword num_points: The number of points on the dispersion curve, equal to the length of the tcp and back_calc arguments. 262 @type num_points: numpy int array of rank [NS][NM][NO] 263 @keyword power: The matrix exponential power array. 264 @type power: numpy int array of rank [NS][NM][NO][ND] 265 """ 266 267 # Once off parameter conversions. 268 dw_AC = dw_AB + dw_BC 269 dwH_AC = dwH_AB + dwH_BC 270 pC = 1.0 - pA - pB 271 pA_pB = pA + pB 272 pA_pC = pA + pC 273 pB_pC = pB + pC 274 k_BA = pA * kex_AB / pA_pB 275 k_AB = pB * kex_AB / pA_pB 276 if pB_pC != 0.0: 277 k_CB = pB * kex_BC / pB_pC 278 k_BC = pC * kex_BC / pB_pC 279 elif pB == 0.0 and pC == 0.0: 280 k_CB = 0.0 281 k_BC = 0.0 282 elif pB == 0.0: 283 k_CB = 0.0 284 k_BC = 1e100 285 elif pC == 0.0: 286 k_CB = 1e100 287 k_BC = 0.0 288 else: 289 k_CB = 1e100 290 k_BC = 1e100 291 k_CA = pA * kex_AC / pA_pC 292 k_AC = pC * kex_AC / pA_pC 293 294 # This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 295 M0[0] = pA 296 M0[1] = pB 297 M0[2] = pC 298 299 # Extract shape of experiment. 300 NS, NM, NO = num_points.shape 301 302 # Populate the m1 and m2 matrices (only once per function call for speed). 303 # D+ matrix component. 304 m1_mat = rmmq_3site_rankN(R20A=R20A, R20B=R20B, R20C=R20C, dw_AB=-dw_AB - dwH_AB, dw_AC=-dw_AC - dwH_AC, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC, k_CB=k_CB, k_AC=k_AC, k_CA=k_CA, tcp=tcp) 305 # Z- matrix component. 306 m2_mat = rmmq_3site_rankN(R20A=R20A, R20B=R20B, R20C=R20C, dw_AB=dw_AB - dwH_AB, dw_AC=dw_AC - dwH_AC, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC, k_CB=k_CB, k_AC=k_AC, k_CA=k_CA, tcp=tcp) 307 308 # The M1 and M2 matrices. 309 # Equivalent to D+. 310 M1_mat = matrix_exponential(m1_mat) 311 # Equivalent to Z-. 312 M2_mat = matrix_exponential(m2_mat) 313 314 # The complex conjugates M1* and M2* 315 # Equivalent to D+*. 316 M1_star_mat = conj(M1_mat) 317 # Equivalent to Z-*. 318 M2_star_mat = conj(M2_mat) 319 320 # Repetitive dot products (minimised for speed). 321 M1_M2_mat = einsum('...ij, ...jk', M1_mat, M2_mat) 322 M2_M1_mat = einsum('...ij, ...jk', M2_mat, M1_mat) 323 M1_M2_M2_M1_mat = einsum('...ij, ...jk', M1_M2_mat, M2_M1_mat) 324 M2_M1_M1_M2_mat = einsum('...ij, ...jk', M2_M1_mat, M1_M2_mat) 325 M1_M2_star_mat = einsum('...ij, ...jk', M1_star_mat, M2_star_mat) 326 M2_M1_star_mat = einsum('...ij, ...jk', M2_star_mat, M1_star_mat) 327 M1_M2_M2_M1_star_mat = einsum('...ij, ...jk', M1_M2_star_mat, M2_M1_star_mat) 328 M2_M1_M1_M2_star_mat = einsum('...ij, ...jk', M2_M1_star_mat, M1_M2_star_mat) 329 330 # Loop over spins. 331 for si in range(NS): 332 # Loop over the spectrometer frequencies. 333 for mi in range(NM): 334 # Loop over offsets: 335 for oi in range(NO): 336 # Extract parameters from array. 337 num_points_i = num_points[si, mi, oi] 338 339 # Loop over the time points, back calculating the R2eff values. 340 for i in range(num_points_i): 341 # Extract data from array. 342 power_i = int(power[si, mi, oi, i]) 343 M1_M2_i = M1_M2_mat[si, mi, oi, i] 344 M1_M2_star_i = M1_M2_star_mat[si, mi, oi, i] 345 M2_M1_i = M2_M1_mat[si, mi, oi, i] 346 M2_M1_star_i = M2_M1_star_mat[si, mi, oi, i] 347 M1_M2_M2_M1_i = M1_M2_M2_M1_mat[si, mi, oi, i] 348 M2_M1_M1_M2_star_i = M2_M1_M1_M2_star_mat[si, mi, oi, i] 349 M2_M1_M1_M2_i = M2_M1_M1_M2_mat[si, mi, oi, i] 350 M1_M2_M2_M1_star_i = M1_M2_M2_M1_star_mat[si, mi, oi, i] 351 352 # Special case of 1 CPMG block - the power is zero. 353 if power_i == 1: 354 # M1.M2. 355 A = M1_M2_i 356 357 # M1*.M2*. 358 B = M1_M2_star_i 359 360 # M2.M1. 361 C = M2_M1_i 362 363 # M2*.M1*. 364 D = M2_M1_star_i 365 366 # Matrices for even number of CPMG blocks. 367 elif power_i % 2 == 0: 368 # The power factor (only calculate once). 369 fact = int(floor(power_i / 2)) 370 371 # (M1.M2.M2.M1)^(n/2). 372 A = matrix_power(M1_M2_M2_M1_i, fact) 373 374 # (M2*.M1*.M1*.M2*)^(n/2). 375 B = matrix_power(M2_M1_M1_M2_star_i, fact) 376 377 # (M2.M1.M1.M2)^(n/2). 378 C = matrix_power(M2_M1_M1_M2_i, fact) 379 380 # (M1*.M2*.M2*.M1*)^(n/2). 381 D = matrix_power(M1_M2_M2_M1_star_i, fact) 382 383 # Matrices for odd number of CPMG blocks. 384 else: 385 # The power factor (only calculate once). 386 fact = int(floor((power_i - 1) / 2)) 387 388 # (M1.M2.M2.M1)^((n-1)/2).M1.M2. 389 A = matrix_power(M1_M2_M2_M1_i, fact) 390 A = dot(A, M1_M2_i) 391 392 # (M1*.M2*.M2*.M1*)^((n-1)/2).M1*.M2*. 393 B = matrix_power(M1_M2_M2_M1_star_i, fact) 394 B = dot(B, M1_M2_star_i) 395 396 # (M2.M1.M1.M2)^((n-1)/2).M2.M1. 397 C = matrix_power(M2_M1_M1_M2_i, fact) 398 C = dot(C, M2_M1_i) 399 400 # (M2*.M1*.M1*.M2*)^((n-1)/2).M2*.M1*. 401 D = matrix_power(M2_M1_M1_M2_star_i, fact) 402 D = dot(D, M2_M1_star_i) 403 404 # The next lines calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential. 405 A_B = dot(A, B) 406 C_D = dot(C, D) 407 Mx = dot(dot(F_vector, (A_B + C_D)), M0) 408 Mx = Mx.real / 2.0 409 if Mx <= 0.0 or isNaN(Mx): 410 back_calc[si, mi, oi, i] = 1e99 411 else: 412 back_calc[si, mi, oi, i]= -inv_tcpmg[si, mi, oi, i] * log(Mx / pA)

413 414

415 -def r2eff_ns_mmq_3site_sq_dq_zq(M0=None, F_vector=array([1, 0, 0], float64), R20A=None, R20B=None, R20C=None, pA=None, pB=None, dw_AB=None, dw_BC=None, dwH_AB=None, dwH_BC=None, kex_AB=None, kex_BC=None, kex_AC=None, inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None):

416 """The 3-site numerical solution to the Bloch-McConnell equation for SQ, ZQ, and DQ data. 417 418 The notation used here comes from: 419 420 - Dmitry M. Korzhnev, Philipp Neudecker, Anthony Mittermaier, Vladislav Yu. Orekhov, and Lewis E. Kay (2005). Multiple-site exchange in proteins studied with a suite of six NMR relaxation dispersion experiments: An application to the folding of a Fyn SH3 domain mutant. J. Am. Chem. Soc., 127, 15602-15611. (doi: http://dx.doi.org/10.1021/ja054550e). 421 422 This function calculates and stores the R2eff values. 423 424 425 @keyword M0: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 426 @type M0: numpy float64, rank-1, 7D array 427 @keyword F_vector: The observable magnitisation vector. This defaults to [1, 0] for X observable magnitisation. 428 @type F_vector: numpy rank-1, 3D float64 array 429 @keyword R20A: The transverse, spin-spin relaxation rate for state A. 430 @type R20A: numpy float array of rank [NS][NM][NO][ND] 431 @keyword R20B: The transverse, spin-spin relaxation rate for state B. 432 @type R20B: numpy float array of rank [NS][NM][NO][ND] 433 @keyword R20C: The transverse, spin-spin relaxation rate for state C. 434 @type R20C: numpy float array of rank [NS][NM][NO][ND] 435 @keyword pA: The population of state A. 436 @type pA: float 437 @keyword pB: The population of state B. 438 @type pB: float 439 @keyword dw_AB: The combined chemical exchange difference between states A and B in rad/s. It should be set to dwH for 1H SQ data, dw for heteronuclear SQ data, dwH-dw for ZQ data, and dwH+dw for DQ data. 440 @type dw_AB: numpy float array of rank [NS][NM][NO][ND] 441 @keyword dw_BC: The combined chemical exchange difference between states B and C in rad/s. It should be set to dwH for 1H SQ data, dw for heteronuclear SQ data, dwH-dw for ZQ data, and dwH+dw for DQ data. 442 @type dw_BC: numpy float array of rank [NS][NM][NO][ND] 443 @keyword dwH_AB: Unused - this is simply to match the r2eff_mmq_3site_mq() function arguments. 444 @type dwH_AB: numpy float array of rank [NS][NM][NO][ND] 445 @keyword dwH_BC: Unused - this is simply to match the r2eff_mmq_3site_mq() function arguments. 446 @type dwH_BC: numpy float array of rank [NS][NM][NO][ND] 447 @keyword kex_AB: The exchange rate between sites A and B for 3-site exchange with kex_AB = k_AB + k_BA (rad.s^-1) 448 @type kex_AB: float 449 @keyword kex_BC: The exchange rate between sites A and C for 3-site exchange with kex_AC = k_AC + k_CA (rad.s^-1) 450 @type kex_BC: float 451 @keyword kex_AC: The exchange rate between sites B and C for 3-site exchange with kex_BC = k_BC + k_CB (rad.s^-1) 452 @type kex_AC: float 453 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds). 454 @type inv_tcpmg: numpy float array of rank [NS][NM][NO][ND] 455 @keyword tcp: The tau_CPMG times (1 / 4.nu1). 456 @type tcp: numpy float array of rank [NS][NM][NO][ND] 457 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies. 458 @type back_calc: numpy float array of rank [NS][NM][NO][ND] 459 @keyword num_points: The number of points on the dispersion curve, equal to the length of the tcp and back_calc arguments. 460 @type num_points: numpy int array of rank [NS][NM][NO] 461 @keyword power: The matrix exponential power array. 462 @type power: numpy int array of rank [NS][NM][NO][ND] 463 """ 464 465 # Once off parameter conversions. 466 dw_AC = dw_AB + dw_BC 467 pC = 1.0 - pA - pB 468 pA_pB = pA + pB 469 pA_pC = pA + pC 470 pB_pC = pB + pC 471 k_BA = pA * kex_AB / pA_pB 472 k_AB = pB * kex_AB / pA_pB 473 if pB_pC != 0.0: 474 k_CB = pB * kex_BC / pB_pC 475 k_BC = pC * kex_BC / pB_pC 476 elif pB == 0.0 and pC == 0.0: 477 k_CB = 0.0 478 k_BC = 0.0 479 elif pB == 0.0: 480 k_CB = 0.0 481 k_BC = 1e100 482 elif pC == 0.0: 483 k_CB = 1e100 484 k_BC = 0.0 485 else: 486 k_CB = 1e100 487 k_BC = 1e100 488 k_CA = pA * kex_AC / pA_pC 489 k_AC = pC * kex_AC / pA_pC 490 491 # This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 492 M0[0] = pA 493 M0[1] = pB 494 M0[2] = pC 495 496 # Extract shape of experiment. 497 NS, NM, NO = num_points.shape 498 499 # Populate the m1 and m2 matrices (only once per function call for speed). 500 # D+ matrix component. 501 m1_mat = rmmq_3site_rankN(R20A=R20A, R20B=R20B, R20C=R20C, dw_AB=dw_AB, dw_AC=dw_AC, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC, k_CB=k_CB, k_AC=k_AC, k_CA=k_CA, tcp=tcp) 502 # Z- matrix component. 503 m2_mat = rmmq_3site_rankN(R20A=R20A, R20B=R20B, R20C=R20C, dw_AB=-dw_AB, dw_AC=-dw_AC, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC, k_CB=k_CB, k_AC=k_AC, k_CA=k_CA, tcp=tcp) 504 505 # The A+/- matrices. 506 A_pos_mat = matrix_exponential(m1_mat) 507 A_neg_mat = matrix_exponential(m2_mat) 508 509 # The evolution for one n. 510 evol_block_mat = einsum('...ij, ...jk', A_neg_mat, A_pos_mat) 511 evol_block_mat = einsum('...ij, ...jk', A_neg_mat, evol_block_mat) 512 evol_block_mat = einsum('...ij, ...jk', A_pos_mat, evol_block_mat) 513 514 # Loop over spins. 515 for si in range(NS): 516 # Loop over the spectrometer frequencies. 517 for mi in range(NM): 518 # Loop over offsets: 519 for oi in range(NO): 520 # Extract parameters from array. 521 num_points_i = num_points[si, mi, oi] 522 523 # Loop over the time points, back calculating the R2eff values. 524 for i in range(num_points_i): 525 # Extract data from array. 526 power_i = int(power[si, mi, oi, i]) 527 evol_block_i = evol_block_mat[si, mi, oi, i] 528 529 # The full evolution. 530 evol = matrix_power(evol_block_i, power_i) 531 532 # The next lines calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential. 533 Mx = dot(F_vector, dot(evol, M0)) 534 Mx = Mx.real 535 if Mx <= 0.0 or isNaN(Mx): 536 back_calc[si, mi, oi, i] = 1e99 537 else: 538 back_calc[si, mi, oi, i] = -inv_tcpmg[si, mi, oi, i] * log(Mx / pA)

539

Source Code for Module lib.dispersion.ns_mmq_3site