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
28
29 """The numerical fit of 2-site Bloch-McConnell equations for CPMG-type experiments, the U{NS CPMG 2-site star<http://wiki.nmr-relax.com/NS_CPMG_2-site_star>} and U{NS CPMG 2-site star full<http://wiki.nmr-relax.com/NS_CPMG_2-site_star_full>} models.
30
31 Description
32 ===========
33
34 The function uses an explicit matrix that contains relaxation, exchange and chemical shift terms. It does the 180deg pulses in the CPMG train with conjugate complex matrices. The approach of Bloch-McConnell can be found in chapter 3.1 of Palmer, A. G. 2004 I{Chem. Rev.}, B{104}, 3623-3640. This function was written, initially in MATLAB, in 2010.
35
36
37 Code origin
38 ===========
39
40 The code was submitted at U{http://thread.gmane.org/gmane.science.nmr.relax.devel/4132} by Paul Schanda.
41
42
43 Links
44 =====
45
46 More information on the NS CPMG 2-site star model can be found in the:
47
48 - U{relax wiki<http://wiki.nmr-relax.com/NS_CPMG_2-site_star>},
49 - U{relax manual<http://www.nmr-relax.com/manual/The_reduced_NS_2_site_star_CPMG_model.html>},
50 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_CPMG_2-site_star>}.
51
52 More information on the NS CPMG 2-site star full model can be found in the:
53
54 - U{relax wiki<http://wiki.nmr-relax.com/NS_CPMG_2-site_star_full>},
55 - U{relax manual<http://www.nmr-relax.com/manual/The_full_NS_2_site_star_CPMG_model.html>},
56 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_CPMG_2-site_star_full>}.
57 """
58
59
60 from numpy import add, array, conj, dot, einsum, fabs, float64, isfinite, log, min, multiply, sum
61 from numpy.ma import fix_invalid, masked_where
62 from numpy.linalg import matrix_power
63
64
65 from lib.float import isNaN
66 from lib.dispersion.matrix_exponential import matrix_exponential
67
68
69 m_r20a = array([
70 [-1, 0],
71 [ 0, 0]], float64)
72
73 m_r20b = array([
74 [ 0, 0],
75 [ 0, -1]], float64)
76
77 m_k_AB = array([
78 [-1, 0],
79 [ 1, 0]], float64)
80
81 m_k_BA = array([
82 [ 0, 1],
83 [ 0, -1]], float64)
84
85 m_dw = array([
86 [ 0, 0],
87 [ 0, 1]], float64)
88
89
90 -def rcpmg_star_rankN(R2A=None, R2B=None, dw=None, k_AB=None, k_BA=None, tcp=None):
91 """Definition of the exchange matrix, for rank [NE][NS][NM][NO][ND][2][2].
92
93 @keyword R2A: The transverse, spin-spin relaxation rate for state A.
94 @type R2A: numpy float array of rank [NE][NS][NM][NO][ND]
95 @keyword R2B: The transverse, spin-spin relaxation rate for state B.
96 @type R2B: numpy float array of rank [NE][NS][NM][NO][ND]
97 @keyword dw: The chemical exchange difference between states A and B in rad/s.
98 @type dw: numpy float array of rank [NE][NS][NM][NO][ND]
99 @keyword k_AB: The forward exchange rate from state A to state B.
100 @type k_AB: float
101 @keyword k_BA: The reverse exchange rate from state B to state A.
102 @type k_BA: float
103 @keyword tcp: The tau_CPMG times (1 / 4.nu1).
104 @type tcp: numpy float array of rank [NE][NS][NM][NO][ND]
105 @return: The relaxation matrix R and complex conjugate cR2.
106 @rtype: numpy float array of rank [NE][NS][NM][NO][ND][2][2]
107 """
108
109
110 r20a_tcp = R2A * tcp
111 r20b_tcp = R2B * tcp
112 k_AB_tcp = k_AB * tcp
113 k_BA_tcp = k_BA * tcp
114
115 dw_tcp_C = dw * tcp * -1j
116
117
118
119
120
121
122
123 m_r20a_tcp = multiply.outer( r20a_tcp, m_r20a )
124 m_r20b_tcp = multiply.outer( r20b_tcp, m_r20b )
125
126
127 Rr_mat = (m_r20a_tcp + m_r20b_tcp)
128
129
130
131
132
133
134
135
136
137 m_k_AB_tcp = multiply.outer( k_AB_tcp, m_k_AB )
138 m_k_BA_tcp = multiply.outer( k_BA_tcp, m_k_BA )
139
140
141 Rex_mat = (m_k_AB_tcp + m_k_BA_tcp)
142
143
144
145
146
147
148 m_dw_tcp_C = multiply.outer( dw_tcp_C, m_dw )
149
150
151 RCS_mat = m_dw_tcp_C
152
153
154 R_mat = add(Rr_mat, Rex_mat)
155 R_mat = add(R_mat, RCS_mat)
156
157
158 cR2_mat = conj(R_mat) * 2.0
159
160
161 return R_mat, cR2_mat, Rr_mat, Rex_mat, RCS_mat
162
163
164 -def r2eff_ns_cpmg_2site_star(M0=None, 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):
165 """The 2-site numerical solution to the Bloch-McConnell equation using complex conjugate matrices.
166
167 This function calculates and stores the R2eff values.
168
169
170 @keyword M0: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations.
171 @type M0: numpy float64, rank-1, 2D array
172 @keyword r20a: The R2 value for state A in the absence of exchange.
173 @type r20a: numpy float array of rank [NE][NS][NM][NO][ND]
174 @keyword r20b: The R2 value for state B in the absence of exchange.
175 @type r20b: numpy float array of rank [NE][NS][NM][NO][ND]
176 @keyword pA: The population of state A.
177 @type pA: float
178 @keyword dw: The chemical exchange difference between states A and B in rad/s.
179 @type dw: numpy float array of rank [NE][NS][NM][NO][ND]
180 @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.
181 @type dw_orig: numpy float array of rank-1
182 @keyword kex: The kex parameter value (the exchange rate in rad/s).
183 @type kex: float
184 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds).
185 @type inv_tcpmg: numpy float array of rank [NE][NS][NM][NO][ND]
186 @keyword tcp: The tau_CPMG times (1 / 4.nu1).
187 @type tcp: numpy float array of rank [NE][NS][NM][NO][ND]
188 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies.
189 @type back_calc: numpy float array of rank [NE][NS][NM][NO][ND]
190 @keyword num_points: The number of points on the dispersion curve, equal to the length of the tcp and back_calc arguments.
191 @type num_points: numpy int array of rank [NE][NS][NM][NO]
192 @keyword power: The matrix exponential power array.
193 @type power: numpy int array of rank [NE][NS][NM][NO][ND]
194 """
195
196
197 t_dw_zero = False
198
199
200 if pA == 1.0 or kex == 0.0:
201 back_calc[:] = r20a
202 return
203
204
205 if min(fabs(dw_orig)) == 0.0:
206 t_dw_zero = True
207 mask_dw_zero = masked_where(dw == 0.0, dw)
208
209
210 pB = 1.0 - pA
211 k_BA = pA * kex
212 k_AB = pB * kex
213
214
215 M0[0] = pA
216 M0[1] = pB
217
218
219 NE, NS, NM, NO, ND = back_calc.shape
220
221
222 R_mat, cR2_mat, Rr_mat, Rex_mat, RCS_mat = rcpmg_star_rankN(R2A=r20a, R2B=r20b, dw=dw, k_AB=k_AB, k_BA=k_BA, tcp=tcp)
223
224
225
226 eR_mat = matrix_exponential(R_mat)
227 ecR2_mat = matrix_exponential(cR2_mat)
228
229
230
231 prop_2_mat = evolution_matrix_mat = einsum('...ij, ...jk', eR_mat, ecR2_mat)
232 prop_2_mat = evolution_matrix_mat = einsum('...ij, ...jk', prop_2_mat, eR_mat)
233
234
235 for si in range(NS):
236
237 for mi in range(NM):
238
239 num_points_si_mi = int(num_points[0, si, mi, 0])
240
241
242 for di in range(num_points_si_mi):
243
244 power_si_mi_di = int(power[0, si, mi, 0, di])
245
246
247 prop_2_i = prop_2_mat[0, si, mi, 0, di]
248
249
250 prop_total = matrix_power(prop_2_i, power_si_mi_di)
251
252
253 Moft = dot(prop_total, M0)
254
255
256 Mx = Moft[0].real / M0[0]
257 if Mx <= 0.0 or isNaN(Mx):
258 back_calc[0, si, mi, 0, di] = 1e99
259 else:
260 back_calc[0, si, mi, 0, di]= -inv_tcpmg[0, si, mi, 0, di] * log(Mx)
261
262
263
264 if t_dw_zero:
265 back_calc[mask_dw_zero.mask] = r20a[mask_dw_zero.mask]
266
267
268
269 if not isfinite(sum(back_calc)):
270
271 fix_invalid(back_calc, copy=False, fill_value=1e100)
272