lib.dispersion.ns_mmq

1 ############################################################################### 2 # # 3 # Copyright (C) 2013 Mathilde Lescanne # 4 # Copyright (C) 2013 Dominique Marion # 5 # Copyright (C) 2013-2014 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 k_CB = pB * kex_BC / pB_pC 277 k_BC = pC * kex_BC / pB_pC 278 k_CA = pA * kex_AC / pA_pC 279 k_AC = pC * kex_AC / pA_pC 280 281 # This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 282 M0[0] = pA 283 M0[1] = pB 284 M0[2] = pC 285 286 # Extract shape of experiment. 287 NS, NM, NO = num_points.shape 288 289 # Populate the m1 and m2 matrices (only once per function call for speed). 290 # D+ matrix component. 291 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) 292 # Z- matrix component. 293 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) 294 295 # The M1 and M2 matrices. 296 # Equivalent to D+. 297 M1_mat = matrix_exponential(m1_mat) 298 # Equivalent to Z-. 299 M2_mat = matrix_exponential(m2_mat) 300 301 # The complex conjugates M1* and M2* 302 # Equivalent to D+*. 303 M1_star_mat = conj(M1_mat) 304 # Equivalent to Z-*. 305 M2_star_mat = conj(M2_mat) 306 307 # Repetitive dot products (minimised for speed). 308 M1_M2_mat = einsum('...ij, ...jk', M1_mat, M2_mat) 309 M2_M1_mat = einsum('...ij, ...jk', M2_mat, M1_mat) 310 M1_M2_M2_M1_mat = einsum('...ij, ...jk', M1_M2_mat, M2_M1_mat) 311 M2_M1_M1_M2_mat = einsum('...ij, ...jk', M2_M1_mat, M1_M2_mat) 312 M1_M2_star_mat = einsum('...ij, ...jk', M1_star_mat, M2_star_mat) 313 M2_M1_star_mat = einsum('...ij, ...jk', M2_star_mat, M1_star_mat) 314 M1_M2_M2_M1_star_mat = einsum('...ij, ...jk', M1_M2_star_mat, M2_M1_star_mat) 315 M2_M1_M1_M2_star_mat = einsum('...ij, ...jk', M2_M1_star_mat, M1_M2_star_mat) 316 317 # Loop over spins. 318 for si in range(NS): 319 # Loop over the spectrometer frequencies. 320 for mi in range(NM): 321 # Loop over offsets: 322 for oi in range(NO): 323 # Extract parameters from array. 324 num_points_i = num_points[si, mi, oi] 325 326 # Loop over the time points, back calculating the R2eff values. 327 for i in range(num_points_i): 328 # Extract data from array. 329 power_i = int(power[si, mi, oi, i]) 330 M1_M2_i = M1_M2_mat[si, mi, oi, i] 331 M1_M2_star_i = M1_M2_star_mat[si, mi, oi, i] 332 M2_M1_i = M2_M1_mat[si, mi, oi, i] 333 M2_M1_star_i = M2_M1_star_mat[si, mi, oi, i] 334 M1_M2_M2_M1_i = M1_M2_M2_M1_mat[si, mi, oi, i] 335 M2_M1_M1_M2_star_i = M2_M1_M1_M2_star_mat[si, mi, oi, i] 336 M2_M1_M1_M2_i = M2_M1_M1_M2_mat[si, mi, oi, i] 337 M1_M2_M2_M1_star_i = M1_M2_M2_M1_star_mat[si, mi, oi, i] 338 339 # Special case of 1 CPMG block - the power is zero. 340 if power_i == 1: 341 # M1.M2. 342 A = M1_M2_i 343 344 # M1*.M2*. 345 B = M1_M2_star_i 346 347 # M2.M1. 348 C = M2_M1_i 349 350 # M2*.M1*. 351 D = M2_M1_star_i 352 353 # Matrices for even number of CPMG blocks. 354 elif power_i % 2 == 0: 355 # The power factor (only calculate once). 356 fact = int(floor(power_i / 2)) 357 358 # (M1.M2.M2.M1)^(n/2). 359 A = matrix_power(M1_M2_M2_M1_i, fact) 360 361 # (M2*.M1*.M1*.M2*)^(n/2). 362 B = matrix_power(M2_M1_M1_M2_star_i, fact) 363 364 # (M2.M1.M1.M2)^(n/2). 365 C = matrix_power(M2_M1_M1_M2_i, fact) 366 367 # (M1*.M2*.M2*.M1*)^(n/2). 368 D = matrix_power(M1_M2_M2_M1_star_i, fact) 369 370 # Matrices for odd number of CPMG blocks. 371 else: 372 # The power factor (only calculate once). 373 fact = int(floor((power_i - 1) / 2)) 374 375 # (M1.M2.M2.M1)^((n-1)/2).M1.M2. 376 A = matrix_power(M1_M2_M2_M1_i, fact) 377 A = dot(A, M1_M2_i) 378 379 # (M1*.M2*.M2*.M1*)^((n-1)/2).M1*.M2*. 380 B = matrix_power(M1_M2_M2_M1_star_i, fact) 381 B = dot(B, M1_M2_star_i) 382 383 # (M2.M1.M1.M2)^((n-1)/2).M2.M1. 384 C = matrix_power(M2_M1_M1_M2_i, fact) 385 C = dot(C, M2_M1_i) 386 387 # (M2*.M1*.M1*.M2*)^((n-1)/2).M2*.M1*. 388 D = matrix_power(M2_M1_M1_M2_star_i, fact) 389 D = dot(D, M2_M1_star_i) 390 391 # The next lines calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential. 392 A_B = dot(A, B) 393 C_D = dot(C, D) 394 Mx = dot(dot(F_vector, (A_B + C_D)), M0) 395 Mx = Mx.real / 2.0 396 if Mx <= 0.0 or isNaN(Mx): 397 back_calc[si, mi, oi, i] = 1e99 398 else: 399 back_calc[si, mi, oi, i]= -inv_tcpmg[si, mi, oi, i] * log(Mx / pA)

400 401

402 -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):

403 """The 3-site numerical solution to the Bloch-McConnell equation for SQ, ZQ, and DQ data. 404 405 The notation used here comes from: 406 407 - 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). 408 409 This function calculates and stores the R2eff values. 410 411 412 @keyword M0: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 413 @type M0: numpy float64, rank-1, 7D array 414 @keyword F_vector: The observable magnitisation vector. This defaults to [1, 0] for X observable magnitisation. 415 @type F_vector: numpy rank-1, 3D float64 array 416 @keyword R20A: The transverse, spin-spin relaxation rate for state A. 417 @type R20A: numpy float array of rank [NS][NM][NO][ND] 418 @keyword R20B: The transverse, spin-spin relaxation rate for state B. 419 @type R20B: numpy float array of rank [NS][NM][NO][ND] 420 @keyword R20C: The transverse, spin-spin relaxation rate for state C. 421 @type R20C: numpy float array of rank [NS][NM][NO][ND] 422 @keyword pA: The population of state A. 423 @type pA: float 424 @keyword pB: The population of state B. 425 @type pB: float 426 @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. 427 @type dw_AB: numpy float array of rank [NS][NM][NO][ND] 428 @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. 429 @type dw_BC: numpy float array of rank [NS][NM][NO][ND] 430 @keyword dwH_AB: Unused - this is simply to match the r2eff_mmq_3site_mq() function arguments. 431 @type dwH_AB: numpy float array of rank [NS][NM][NO][ND] 432 @keyword dwH_BC: Unused - this is simply to match the r2eff_mmq_3site_mq() function arguments. 433 @type dwH_BC: numpy float array of rank [NS][NM][NO][ND] 434 @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) 435 @type kex_AB: float 436 @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) 437 @type kex_BC: float 438 @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) 439 @type kex_AC: float 440 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds). 441 @type inv_tcpmg: numpy float array of rank [NS][NM][NO][ND] 442 @keyword tcp: The tau_CPMG times (1 / 4.nu1). 443 @type tcp: numpy float array of rank [NS][NM][NO][ND] 444 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies. 445 @type back_calc: numpy float array of rank [NS][NM][NO][ND] 446 @keyword num_points: The number of points on the dispersion curve, equal to the length of the tcp and back_calc arguments. 447 @type num_points: numpy int array of rank [NS][NM][NO] 448 @keyword power: The matrix exponential power array. 449 @type power: numpy int array of rank [NS][NM][NO][ND] 450 """ 451 452 # Once off parameter conversions. 453 dw_AC = dw_AB + dw_BC 454 pC = 1.0 - pA - pB 455 pA_pB = pA + pB 456 pA_pC = pA + pC 457 pB_pC = pB + pC 458 k_BA = pA * kex_AB / pA_pB 459 k_AB = pB * kex_AB / pA_pB 460 k_CB = pB * kex_BC / pB_pC 461 k_BC = pC * kex_BC / pB_pC 462 k_CA = pA * kex_AC / pA_pC 463 k_AC = pC * kex_AC / pA_pC 464 465 # This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 466 M0[0] = pA 467 M0[1] = pB 468 M0[2] = pC 469 470 # Extract shape of experiment. 471 NS, NM, NO = num_points.shape 472 473 # Populate the m1 and m2 matrices (only once per function call for speed). 474 # D+ matrix component. 475 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) 476 # Z- matrix component. 477 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) 478 479 # The A+/- matrices. 480 A_pos_mat = matrix_exponential(m1_mat) 481 A_neg_mat = matrix_exponential(m2_mat) 482 483 # The evolution for one n. 484 evol_block_mat = einsum('...ij, ...jk', A_neg_mat, A_pos_mat) 485 evol_block_mat = einsum('...ij, ...jk', A_neg_mat, evol_block_mat) 486 evol_block_mat = einsum('...ij, ...jk', A_pos_mat, evol_block_mat) 487 488 # Loop over spins. 489 for si in range(NS): 490 # Loop over the spectrometer frequencies. 491 for mi in range(NM): 492 # Loop over offsets: 493 for oi in range(NO): 494 # Extract parameters from array. 495 num_points_i = num_points[si, mi, oi] 496 497 # Loop over the time points, back calculating the R2eff values. 498 for i in range(num_points_i): 499 # Extract data from array. 500 power_i = int(power[si, mi, oi, i]) 501 evol_block_i = evol_block_mat[si, mi, oi, i] 502 503 # The full evolution. 504 evol = matrix_power(evol_block_i, power_i) 505 506 # The next lines calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential. 507 Mx = dot(F_vector, dot(evol, M0)) 508 Mx = Mx.real 509 if Mx <= 0.0 or isNaN(Mx): 510 back_calc[si, mi, oi, i] = 1e99 511 else: 512 back_calc[si, mi, oi, i] = -inv_tcpmg[si, mi, oi, i] * log(Mx / pA)

513

Source Code for Module lib.dispersion.ns_mmq_3site