1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27 """The numerical fit of 2-site Bloch-McConnell equations for CPMG-type experiments, the U{NS CPMG 2-site 3D<http://wiki.nmr-relax.com/NS_CPMG_2-site_3D>} and U{NS CPMG 2-site 3D full<http://wiki.nmr-relax.com/NS_CPMG_2-site_3D_full>} models.
28
29 Description
30 ===========
31
32 The function uses an explicit matrix that contains relaxation, exchange and chemical shift terms. It does the 180deg pulses in the CPMG train. The approach of Bloch-McConnell can be found in chapter 3.1 of Palmer, A. G. Chem Rev 2004, 104, 3623-3640. This function was written, initially in MATLAB, in 2010.
33
34
35 Code origin
36 ===========
37
38 This is the model of the numerical solution for the 2-site Bloch-McConnell equations. It originates as optimization function number 1 from the fitting_main_kex.py script from Mathilde Lescanne, Paul Schanda, and Dominique Marion (see U{http://thread.gmane.org/gmane.science.nmr.relax.devel/4138}, U{https://web.archive.org/web/https://gna.org/task/?7712#comment2} and U{https://web.archive.org/web/https://gna.org/support/download.php?file_id=18262}).
39
40
41 Links
42 =====
43
44 More information on the NS CPMG 2-site 3D model can be found in the:
45
46 - U{relax wiki<http://wiki.nmr-relax.com/NS_CPMG_2-site_3D>},
47 - U{relax manual<http://www.nmr-relax.com/manual/The_reduced_NS_2_site_3D_CPMG_model.html>},
48 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_CPMG_2-site_3D>}.
49
50 More information on the NS CPMG 2-site 3D full model can be found in the:
51
52 - U{relax wiki<http://wiki.nmr-relax.com/NS_CPMG_2-site_3D_full>},
53 - U{relax manual<http://www.nmr-relax.com/manual/The_full_NS_2_site_3D_CPMG_model.html>},
54 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_CPMG_2-site_3D_full>}.
55 """
56
57
58 from numpy import array, dot, fabs, float64, einsum, isfinite, log, min, multiply, rollaxis, sum
59 from numpy.linalg import matrix_power
60 from numpy.ma import fix_invalid, masked_where
61
62
63 from lib.float import isNaN
64 from lib.dispersion.matrix_exponential import matrix_exponential
65
66
67 m_r10a = array([
68 [0, 0, 0, 0, 0, 0, 0],
69 [0, 0, 0, 0, 0, 0, 0],
70 [0, 0, 0, 0, 0, 0, 0],
71 [1, 0, 0, -1, 0, 0, 0],
72 [0, 0, 0, 0, 0, 0, 0],
73 [0, 0, 0, 0, 0, 0, 0],
74 [0, 0, 0, 0, 0, 0, 0]], float64)
75
76 m_pA = array([
77 [0, 0, 0, 0, 0, 0, 0],
78 [0, 0, 0, 0, 0, 0, 0],
79 [0, 0, 0, 0, 0, 0, 0],
80 [2, 0, 0, 0, 0, 0, 0],
81 [0, 0, 0, 0, 0, 0, 0],
82 [0, 0, 0, 0, 0, 0, 0],
83 [0, 0, 0, 0, 0, 0, 0]], float64)
84
85 m_r10b = array([
86 [0, 0, 0, 0, 0, 0, 0],
87 [0, 0, 0, 0, 0, 0, 0],
88 [0, 0, 0, 0, 0, 0, 0],
89 [0, 0, 0, 0, 0, 0, 0],
90 [0, 0, 0, 0, 0, 0, 0],
91 [0, 0, 0, 0, 0, 0, 0],
92 [1, 0, 0, 0, 0, 0, -1]], float64)
93
94 m_pB = array([
95 [0, 0, 0, 0, 0, 0, 0],
96 [0, 0, 0, 0, 0, 0, 0],
97 [0, 0, 0, 0, 0, 0, 0],
98 [0, 0, 0, 0, 0, 0, 0],
99 [0, 0, 0, 0, 0, 0, 0],
100 [0, 0, 0, 0, 0, 0, 0],
101 [2, 0, 0, 0, 0, 0, 0]], float64)
102
103 m_r20a = array([
104 [0, 0, 0, 0, 0, 0, 0],
105 [0, -1, 0, 0, 0, 0, 0],
106 [0, 0, -1, 0, 0, 0, 0],
107 [0, 0, 0, 0, 0, 0, 0],
108 [0, 0, 0, 0, 0, 0, 0],
109 [0, 0, 0, 0, 0, 0, 0],
110 [0, 0, 0, 0, 0, 0, 0]], float64)
111
112 m_r20b = array([
113 [0, 0, 0, 0, 0, 0, 0],
114 [0, 0, 0, 0, 0, 0, 0],
115 [0, 0, 0, 0, 0, 0, 0],
116 [0, 0, 0, 0, 0, 0, 0],
117 [0, 0, 0, 0, -1, 0, 0],
118 [0, 0, 0, 0, 0, -1, 0],
119 [0, 0, 0, 0, 0, 0, 0]], float64)
120
121 m_k_AB = array([
122 [0, 0, 0, 0, 0, 0, 0],
123 [0, -1, 0, 0, 0, 0, 0],
124 [0, 0, -1, 0, 0, 0, 0],
125 [0, 0, 0, -1, 0, 0, 0],
126 [0, 1, 0, 0, 0, 0, 0],
127 [0, 0, 1, 0, 0, 0, 0],
128 [0, 0, 0, 1, 0, 0, 0]], float64)
129
130 m_k_BA = array([
131 [0, 0, 0, 0, 0, 0, 0],
132 [0, 0, 0, 0, 1, 0, 0],
133 [0, 0, 0, 0, 0, 1, 0],
134 [0, 0, 0, 0, 0, 0, 1],
135 [0, 0, 0, 0, -1, 0, 0],
136 [0, 0, 0, 0, 0, -1, 0],
137 [0, 0, 0, 0, 0, 0, -1]], float64)
138
139 m_wA = array([
140 [0, 0, 0, 0, 0, 0, 0],
141 [0, 0, -1, 0, 0, 0, 0],
142 [0, 1, 0, 0, 0, 0, 0],
143 [0, 0, 0, 0, 0, 0, 0],
144 [0, 0, 0, 0, 0, 0, 0],
145 [0, 0, 0, 0, 0, 0, 0],
146 [0, 0, 0, 0, 0, 0, 0]], float64)
147
148 m_wB = array([
149 [0, 0, 0, 0, 0, 0, 0],
150 [0, 0, 0, 0, 0, 0, 0],
151 [0, 0, 0, 0, 0, 0, 0],
152 [0, 0, 0, 0, 0, 0, 0],
153 [0, 0, 0, 0, 0, -1, 0],
154 [0, 0, 0, 0, 1, 0, 0],
155 [0, 0, 0, 0, 0, 0, 0]], float64)
156
157
158 -def rcpmg_3d_rankN(R1A=None, R1B=None, R2A=None, R2B=None, pA=None, pB=None, dw=None, k_AB=None, k_BA=None, tcp=None):
159 """Definition of the 3D exchange matrix, for rank [NE][NS][NM][NO][ND][7][7].
160
161 @keyword R1A: The longitudinal, spin-lattice relaxation rate for state A.
162 @type R1A: numpy float array of rank [NE][NS][NM][NO][ND]
163 @keyword R1B: The longitudinal, spin-lattice relaxation rate for state B.
164 @type R1B: numpy float array of rank [NE][NS][NM][NO][ND]
165 @keyword R2A: The transverse, spin-spin relaxation rate for state A.
166 @type R2A: numpy float array of rank [NE][NS][NM][NO][ND]
167 @keyword R2B: The transverse, spin-spin relaxation rate for state B.
168 @type R2B: numpy float array of rank [NE][NS][NM][NO][ND]
169 @keyword pA: The population of state A.
170 @type pA: float
171 @keyword pB: The population of state B.
172 @type pB: float
173 @keyword dw: The chemical exchange difference between states A and B in rad/s.
174 @type dw: numpy float array of rank [NE][NS][NM][NO][ND]
175 @keyword k_AB: The forward exchange rate from state A to state B.
176 @type k_AB: float
177 @keyword k_BA: The reverse exchange rate from state B to state A.
178 @type k_BA: float
179 @keyword tcp: The tau_CPMG times (1 / 4.nu1).
180 @type tcp: numpy float array of rank [NE][NS][NM][NO][ND]
181 @return: The relaxation matrix.
182 @rtype: numpy float array of rank [NE][NS][NM][NO][ND][7][7]
183 """
184
185
186 wA = 0.0
187 wB = dw
188
189 r10a_tcp = R1A * tcp
190 r10b_tcp = R1B * tcp
191 r20a_tcp = R2A * tcp
192 r20b_tcp = R2B * tcp
193 pA_tcp = pA * tcp
194 pB_tcp = pB * tcp
195 dw_tcp = dw * tcp
196 k_AB_tcp = k_AB * tcp
197 k_BA_tcp = k_BA * tcp
198 wA_tcp = wA * tcp
199 wB_tcp = wB * tcp
200
201
202
203 m_r10a_tcp = multiply.outer( r10a_tcp, m_r10a )
204 m_pA_tcp = multiply.outer( pA_tcp, m_pA )
205
206 m_r10b_tcp = multiply.outer( r10b_tcp, m_r10b )
207 m_pB_tcp = multiply.outer( pB_tcp, m_pB )
208
209 m_r20a_tcp = multiply.outer( r20a_tcp, m_r20a )
210 m_r20b_tcp = multiply.outer( r20b_tcp, m_r20b )
211
212 m_k_AB_tcp = multiply.outer( k_AB_tcp, m_k_AB )
213 m_k_BA_tcp = multiply.outer( k_BA_tcp, m_k_BA )
214
215 m_wA_tcp = multiply.outer( wA_tcp, m_wA )
216 m_wB_tcp = multiply.outer( wB_tcp, m_wB )
217
218
219 c_mat = (m_r10a_tcp * m_pA_tcp + m_r10b_tcp * m_pB_tcp
220 + m_r20a_tcp + m_r20b_tcp
221 + m_k_AB_tcp + m_k_BA_tcp
222 + m_wA_tcp + m_wB_tcp )
223
224
225 return c_mat
226
227
228 -def r2eff_ns_cpmg_2site_3D(r180x=None, M0=None, M0_T=None, r10a=0.0, r10b=0.0, r20a=None, r20b=None, pA=None, dw=None, dw_orig=None, kex=None, inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None):
229 """The 2-site numerical solution to the Bloch-McConnell equation.
230
231 This function calculates and stores the R2eff values.
232
233
234 @keyword r180x: The X-axis pi-pulse propagator.
235 @type r180x: numpy float64, rank-2, 7D array
236 @keyword M0: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations.
237 @type M0: numpy float array of rank [NE][NS][NM][NO][ND][7][1]
238 @keyword M0_T: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations, where the outer two axis has been swapped for efficient dot operations.
239 @type M0_T: numpy float array of rank [NE][NS][NM][NO][ND][1][7]
240 @keyword r10a: The R1 value for state A.
241 @type r10a: float
242 @keyword r10b: The R1 value for state B.
243 @type r10b: float
244 @keyword r20a: The R2 value for state A in the absence of exchange.
245 @type r20a: numpy float array of rank [NE][NS][NM][NO][ND]
246 @keyword r20b: The R2 value for state B in the absence of exchange.
247 @type r20b: numpy float array of rank [NE][NS][NM][NO][ND]
248 @keyword pA: The population of state A.
249 @type pA: float
250 @keyword dw: The chemical exchange difference between states A and B in rad/s.
251 @type dw: numpy float array of rank [NE][NS][NM][NO][ND]
252 @keyword dw_orig: The chemical exchange difference between states A and B in ppm. This is only for faster checking of zero value, which result in no exchange.
253 @type dw_orig: numpy float array of rank-1
254 @keyword kex: The kex parameter value (the exchange rate in rad/s).
255 @type kex: float
256 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds).
257 @type inv_tcpmg: numpy float array of rank [NE][NS][NM][NO][ND]
258 @keyword tcp: The tau_CPMG times (1 / 4.nu1).
259 @type tcp: numpy float array of rank [NE][NS][NM][NO][ND]
260 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies.
261 @type back_calc: numpy float array of rank [NE][NS][NM][NO][ND]
262 @keyword num_points: The number of points on the dispersion curve, equal to the length of the tcp and back_calc arguments.
263 @type num_points: numpy int array of rank [NE][NS][NM][NO]
264 @keyword power: The matrix exponential power array.
265 @type power: numpy int array of rank [NE][NS][NM][NO][ND]
266 """
267
268
269 t_dw_zero = False
270
271
272 if pA == 1.0 or kex == 0.0:
273 back_calc[:] = r20a
274 return
275
276
277 if min(fabs(dw_orig)) == 0.0:
278 t_dw_zero = True
279 mask_dw_zero = masked_where(dw == 0.0, dw)
280
281
282 pB = 1.0 - pA
283 k_BA = pA * kex
284 k_AB = pB * kex
285
286
287 M0_T[:, :, :, :, :, 0, 1] = pA
288 M0_T[:, :, :, :, :, 0, 4] = pB
289 M0[:, :, :, :, :, 1, 0] = pA
290 M0[:, :, :, :, :, 4, 0] = pB
291
292
293 NE, NS, NM, NO, ND = back_calc.shape
294
295
296 R_mat = rcpmg_3d_rankN(R1A=r10a, R1B=r10b, R2A=r20a, R2B=r20b, pA=pA, pB=pB, dw=dw, k_AB=k_AB, k_BA=k_BA, tcp=tcp)
297
298
299 Rexpo_mat = matrix_exponential(R_mat)
300
301
302
303
304
305 evolution_matrix_mat = einsum('...ij,...jk', Rexpo_mat, r180x)
306 evolution_matrix_mat = einsum('...ij,...jk', evolution_matrix_mat, Rexpo_mat)
307 evolution_matrix_mat = einsum('...ij,...jk', evolution_matrix_mat, evolution_matrix_mat)
308
309
310 evolution_matrix_T_mat = rollaxis(evolution_matrix_mat, 6, 5)
311
312
313 evolution_matrix_T_M0_mat = einsum('...ij,...jk', M0_T, evolution_matrix_T_mat)
314
315
316 for si in range(NS):
317
318 for mi in range(NM):
319
320 num_points_si_mi = int(num_points[0, si, mi, 0])
321
322
323 for di in range(num_points_si_mi):
324
325 inv_tcpmg_si_mi_di = inv_tcpmg[0, si, mi, 0, di]
326 power_si_mi_di = int(power[0, si, mi, 0, di])
327 r20a_si_mi_di = r20a[0, si, mi, 0, di]
328
329
330 Mint_T_i = evolution_matrix_T_M0_mat[0, si, mi, 0, di]
331
332
333 evolution_matrix_T_i = evolution_matrix_T_mat[0, si, mi, 0, di]
334
335
336 l = int(power_si_mi_di-1)
337
338
339 evolution_matrix_T_power_i = matrix_power(evolution_matrix_T_i, l)
340
341
342 Mint_T_i = dot(Mint_T_i, evolution_matrix_T_power_i)
343
344
345 Mx = Mint_T_i[0][1] / pA
346 if Mx <= 0.0 or isNaN(Mx):
347 back_calc[0, si, mi, 0, di] = r20a_si_mi_di
348 else:
349 back_calc[0, si, mi, 0, di] = - inv_tcpmg_si_mi_di * log(Mx)
350
351
352
353 if t_dw_zero:
354 back_calc[mask_dw_zero.mask] = r20a[mask_dw_zero.mask]
355
356
357
358 if not isfinite(sum(back_calc)):
359
360 fix_invalid(back_calc, copy=False, fill_value=1e100)
361