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 expanded<http://wiki.nmr-relax.com/NS_CPMG_2-site_expanded>} model.
30
31 Description
32 ===========
33
34 This function is exact, just as the explicit Bloch-McConnell numerical treatments. It comes from a Maple derivation based on the Bloch-McConnell equations. It is much faster than the numerical Bloch-McConnell solution. It was derived by Nikolai Skrynnikov and is provided with his permission.
35
36
37 Code origin
38 ===========
39
40 The code originates as optimization function number 5 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}).
41
42 Links to the copyright licensing agreements from all authors are:
43
44 - Nikolai Skrynnikov, U{http://article.gmane.org/gmane.science.nmr.relax.devel/4279},
45 - Martin Tollinger, U{http://article.gmane.org/gmane.science.nmr.relax.devel/4276},
46 - Paul Schanda, U{http://article.gmane.org/gmane.science.nmr.relax.devel/4271},
47 - Mathilde Lescanne, U{http://article.gmane.org/gmane.science.nmr.relax.devel/4138},
48 - Dominique Marion, U{http://article.gmane.org/gmane.science.nmr.relax.devel/4157}.
49
50
51 Code evolution
52 --------------
53
54 The complex path of the code from the original Maple to relax can be described as:
55
56 - p3.analytical (Maple input text file at U{https://web.archive.org/web/https://gna.org/task/?7712#comment8}),
57 - Automatically generated FORTRAN,
58 - Manually converted to Matlab by Nikolai (sim_all.tar at U{https://web.archive.org/web/https://gna.org/task/?7712#comment5})
59 - Manually converted to Python by Paul, Mathilde, and Dominique (fitting_main.py at U{https://web.archive.org/web/https://gna.org/task/?7712#comment1})
60 - Converted into Python code for relax (here).
61
62
63 Maple p3.analytical script
64 --------------------------
65
66 For reference, the original Maple script written by Nikolai for the expansion of the equations is::
67
68 with(linalg):
69 with(tensor):
70 #Ka:=30;
71 #Kb:=1200;
72 #dW:=300;
73 #N:=2;
74 #tcp:=0.040/N;
75
76 Ksym:=sqrt(Ka*Kb);
77 #dX:=(Ka-Kb+I*dw)/2; # Ra=Rb
78 dX:=((Ra-Rb)+(Ka-Kb)+I*dw)/2;
79
80 L:=([[-dX, Ksym], [Ksym, dX]]);
81
82 # in the end everything is multiplied by exp(-0.5*(Ra+Rb+Ka+Kb)*(Tc+2*tpalmer))
83 # where 0.5*(Ra+Rb) is the same as Rinf, and (Ka+Kb) is kex.
84
85 y:=eigenvects(L);
86 TP1:=array([[y[1][3][1][1],y[2][3][1][1]],[y[1][3][1][2],y[2][3][1][2]]]);
87 iTP1:=inverse(TP1);
88 P1:=array([[exp(y[1][1]*tcp/2),0],[0,exp(y[2][1]*tcp/2)]]);
89
90 P1palmer:=array([[exp(y[1][1]*tpalmer),0],[0,exp(y[2][1]*tpalmer)]]);
91
92 TP2:=map(z->conj(z),TP1);
93 iTP2:=map(z->conj(z),iTP1);
94 P2:=array([[exp(conj(y[1][1])*tcp),0],[0,exp(conj(y[2][1])*tcp)]]);
95
96 P2palmer:=array([[exp(conj(y[1][1])*tpalmer),0],[0,exp(conj(y[2][1])*tpalmer)]]);
97
98 cP1:=evalm(TP1&*P1&*iTP1);
99 cP2:=evalm(TP2&*P2&*iTP2);
100
101 cP1palmer:=evalm(TP1&*P1palmer&*iTP1);
102 cP2palmer:=evalm(TP2&*P2palmer&*iTP2);
103
104 Ps:=evalm(cP1&*cP2&*cP1);
105 # Ps is symmetric; cf. simplify(Ps[1,2]-Ps[2,1]);
106 Pspalmer:=evalm(cP2palmer&*cP1palmer);
107
108
109 dummy:=array([[a,b],[b,c]]);
110 x:=eigenvects(dummy);
111 TG1:=array([[x[1][3][1][1],x[2][3][1][1]],[x[1][3][1][2],x[2][3][1][2]]]);
112 iTG1:=inverse(TG1);
113 G1:=array([[x[1][1]^(N/4),0],[0,x[2][1]^(N/4)]]);
114 GG1:=evalm(TG1&*G1&*iTG1);
115 GG2:=map(z->conj(z),GG1);
116
117 cGG:=evalm(GG2&*Pspalmer&*GG1);
118
119 #s0:=array([Kb, Ka]);
120 s0:=array([sqrt(Kb),sqrt(Ka)]); # accounts for exchange symmetrization
121 st:=evalm(cGG&*s0);
122 #obs:=(1/(Ka+Kb))*st[1];
123 obs:=(sqrt(Kb)/(Ka+Kb))*st[1]; # accounts for exchange symmetrization
124
125 obs1:=eval(obs,[a=Ps[1,1],b=Ps[1,2],c=Ps[2,2]]);
126 #obs2:=simplify(obs1):
127
128 print(obs1):
129
130 cGGref:=evalm(Pspalmer);
131 stref:=evalm(cGGref&*s0);
132 obsref:=(sqrt(Kb)/(Ka+Kb))*stref[1]; # accounts for exchange symmetrization
133
134 print(obsref):
135
136 writeto(result_test):
137
138 fortran([intensity=obs1, intensity_ref=obsref], optimized):
139
140
141 Matlab sim_all.tar funNikolai.m script
142 --------------------------------------
143
144 Also for reference, the Matlab code from Nikolai and Martin manually converted from the automatically generated FORTRAN from the previous script into the funNikolai.m file is::
145
146 function residual = funNikolai(optpar)
147
148 % extended Carver's equation derived via Maple, Ra-Rb = 0 by Skrynnikov
149
150 global nu_0 x y Rcalc rms nfields
151 global Tc
152
153 Rcalc=zeros(nfields,size(x,2));
154
155 tau_ex=optpar(1);
156 pb=optpar(2);
157
158 pa=1-pb;
159 kex=1/tau_ex;
160 Ka=kex*pb;
161 Kb=kex*pa;
162
163 nu_cpmg=x;
164 tcp=1./(2*nu_cpmg);
165 N=round(Tc./tcp);
166
167 for k=1:nfields
168 dw=2*pi*nu_0(k)*optpar(3);
169 Rinf=optpar(3+k);
170
171 t3 = i;
172 t4 = t3*dw;
173 t5 = Kb^2;
174 t8 = 2*t3*Kb*dw;
175 t10 = 2*Kb*Ka;
176 t11 = dw^2;
177 t14 = 2*t3*Ka*dw;
178 t15 = Ka^2;
179 t17 = sqrt(t5-t8+t10-t11+t14+t15);
180 t21 = exp((-Kb+t4-Ka+t17)*tcp/4);
181 t22 = 1/t17;
182 t28 = exp((-Kb+t4-Ka-t17)*tcp/4);
183 t31 = t21*t22*Ka-t28*t22*Ka;
184 t33 = sqrt(t5+t8+t10-t11-t14+t15);
185 t34 = Kb+t4-Ka+t33;
186 t37 = exp((-Kb-t4-Ka+t33)*tcp/2);
187 t39 = 1/t33;
188 t41 = Kb+t4-Ka-t33;
189 t44 = exp((-Kb-t4-Ka-t33)*tcp/2);
190 t47 = t34*t37*t39/2-t41*t44*t39/2;
191 t49 = Kb-t4-Ka-t17;
192 t51 = t21*t49*t22;
193 t52 = Kb-t4-Ka+t17;
194 t54 = t28*t52*t22;
195 t55 = -t51+t54;
196 t60 = t37*t39*Ka-t44*t39*Ka;
197 t62 = t31.*t47+t55.*t60/2;
198 t63 = 1/Ka;
199 t68 = -t52*t63*t51/2+t49*t63*t54/2;
200 t69 = t62.*t68/2;
201 t72 = t37*t41*t39;
202 t76 = t44*t34*t39;
203 t78 = -t34*t63*t72/2+t41*t63*t76/2;
204 t80 = -t72+t76;
205 t82 = t31.*t78/2+t55.*t80/4;
206 t83 = t82.*t55/2;
207 t88 = t52*t21*t22/2-t49*t28*t22/2;
208 t91 = t88.*t47+t68.*t60/2;
209 t92 = t91.*t88;
210 t95 = t88.*t78/2+t68.*t80/4;
211 t96 = t95.*t31;
212 t97 = t69+t83;
213 t98 = t97.^2;
214 t99 = t92+t96;
215 t102 = t99.^2;
216 t108 = t62.*t88+t82.*t31;
217 t112 = sqrt(t98-2*t99.*t97+t102+4*(t91.*t68/2+t95.*t55/2).*t108);
218 t113 = t69+t83-t92-t96-t112;
219 t115 = N/2;
220 t116 = (t69/2+t83/2+t92/2+t96/2+t112/2).^t115;
221 t118 = 1./t112;
222 t120 = t69+t83-t92-t96+t112;
223 t122 = (t69/2+t83/2+t92/2+t96/2-t112/2).^t115;
224 t127 = 1./t108;
225 t139 = 1/(Ka+Kb)*((-t113.*t116.*t118/2+t120.*t122.*t118/2)*Kb+(-t113.*t127.*t116.*t120.*t118/2+t120.*t127.*t122.*t113.*t118/2)*Ka/2);
226
227 intensity0 = pa; % pA
228 intensity = real(t139)*exp(-Tc*Rinf); % that's "homogeneous" relaxation
229 Rcalc(k,:)=(1/Tc)*log(intensity0./intensity);
230
231 end
232
233 if (size(Rcalc)==size(y))
234 residual=sum(sum((y-Rcalc).^2));
235 rms=sqrt(residual/(size(y,1)*size(y,2)));
236 end
237
238
239 Links
240 =====
241
242 More information on the NS CPMG 2-site expanded model can be found in the:
243
244 - U{relax wiki<http://wiki.nmr-relax.com/NS_CPMG_2-site_expanded>},
245 - U{relax manual<http://www.nmr-relax.com/manual/The_NS_2_site_expanded_CPMG_model.html>},
246 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_CPMG_2-site_expanded>}.
247 """
248
249
250 from numpy import any, exp, isfinite, fabs, power, log, min, sqrt, sum
251 from numpy.ma import fix_invalid, masked_where
252
253
254 -def r2eff_ns_cpmg_2site_expanded(r20=None, pA=None, dw=None, dw_orig=None, kex=None, relax_time=None, inv_relax_time=None, tcp=None, back_calc=None, num_cpmg=None):
255 """The 2-site numerical solution to the Bloch-McConnell equation using complex conjugate matrices.
256
257 This function calculates and stores the R2eff values.
258
259
260 @keyword r20: The R2 value for both states A and B in the absence of exchange.
261 @type r20: numpy float array of rank [NE][NS][NM][NO][ND]
262 @keyword pA: The population of state A.
263 @type pA: float
264 @keyword dw: The chemical exchange difference between states A and B in rad/s.
265 @type dw: numpy float array of rank [NE][NS][NM][NO][ND]
266 @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.
267 @type dw_orig: numpy float array of rank-1
268 @keyword kex: The kex parameter value (the exchange rate in rad/s).
269 @type kex: float
270 @keyword relax_time: The total relaxation time period (in seconds).
271 @type relax_time: numpy float array of rank [NE][NS][NM][NO][ND]
272 @keyword inv_relax_time: The inverse of the total relaxation time period (in inverse seconds).
273 @type inv_relax_time: numpy float array of rank [NE][NS][NM][NO][ND]
274 @keyword tcp: The tau_CPMG times (1 / 4.nu1).
275 @type tcp: numpy float array of rank [NE][NS][NM][NO][ND]
276 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies.
277 @type back_calc: numpy float array of rank [NE][NS][NM][NO][ND]
278 @keyword num_cpmg: The array of numbers of CPMG blocks.
279 @type num_cpmg: numpy int16 array of rank [NE][NS][NM][NO][ND]
280 """
281
282
283 t_dw_zero = False
284 t_t108_zero = False
285 t_t112_zero = False
286
287
288 if pA == 1.0 or kex == 0.0:
289 back_calc[:] = r20
290 return
291
292
293 if min(fabs(dw_orig)) == 0.0:
294 t_dw_zero = True
295 mask_dw_zero = masked_where(dw == 0.0, dw)
296
297
298 pB = 1.0 - pA
299 k_BA = pA * kex
300 k_AB = pB * kex
301
302
303 half_tcp = 0.5 * tcp
304 k_AB_plus_k_BA = k_AB + k_BA
305 k_BA_minus_k_AB = k_BA - k_AB
306
307
308 t3 = 1.j
309 t4 = t3 * dw
310 two_t4 = 2.0 * t4
311 t5 = k_BA**2
312 t8 = two_t4 * k_BA
313 t10 = 2.0 * k_BA * k_AB
314 t11 = dw**2
315 t14 = two_t4 * k_AB
316 t15 = k_AB**2
317 t5_t10_t11_t15 = t5 + t10 - t11 + t15
318 t8_t14 = t8 - t14
319 t17 = sqrt(t5_t10_t11_t15 - t8_t14)
320
321 k_AB_plus_k_BA_minus_t4 = k_AB_plus_k_BA - t4
322 t21 = exp((t17 - k_AB_plus_k_BA_minus_t4) * half_tcp)
323 t22 = 1.0/t17
324 t28 = exp(-(t17 + k_AB_plus_k_BA_minus_t4) * half_tcp)
325 t31 = t22*k_AB * (t21 - t28)
326 t33 = sqrt(t5_t10_t11_t15 + t8_t14)
327
328 k_AB_plus_k_BA_plus_t4 = k_AB_plus_k_BA + t4
329 k_BA_minus_k_AB_plus_t4 = k_BA_minus_k_AB + t4
330 t34 = k_BA_minus_k_AB_plus_t4 + t33
331 t37 = exp((t33 - k_AB_plus_k_BA_plus_t4) * tcp)
332 t39 = 1.0/t33
333 t41 = k_BA_minus_k_AB_plus_t4 - t33
334 t44 = exp(-(t33 + k_AB_plus_k_BA_plus_t4) * tcp)
335 t47 = 0.5*t39 * (t34*t37 - t41*t44)
336
337 k_BA_minus_k_AB_minus_t4 = k_BA_minus_k_AB - t4
338 t49 = k_BA_minus_k_AB_minus_t4 - t17
339 t51 = t21 * t49 * t22
340 t52 = k_BA_minus_k_AB_minus_t4 + t17
341 t54 = t28 * t52 * t22
342 t55 = -t51 + t54
343 t60 = 0.5*t39*k_AB * (t37 - t44)
344 t62 = t31*t47 + t55*t60
345 t63 = 1.0/k_AB
346 t68 = 0.5*t63 * (t49*t54 - t52*t51)
347 t69 = 0.5*t62 * t68
348 t72 = t37 * t41 * t39
349 t76 = t44 * t34 * t39
350 t78 = 0.5*t63 * (t41*t76 - t34*t72)
351 t80 = 0.5 * (t76 - t72)
352 t82 = 0.5 * (t31*t78 + t55*t80)
353 t83 = t82 * t55/2.0
354 t88 = 0.5 * t22 * (t52*t21 - t49*t28)
355 t91 = t88 * t47 + t68*t60
356 t92 = t91 * t88
357 t95 = 0.5 * (t88*t78 + t68*t80)
358 t96 = t95 * t31
359 t97 = t69 + t83
360 t98 = t97**2
361 t99 = t92 + t96
362 t102 = t99**2
363 t108 = t62 * t88 + t82 * t31
364
365
366 mask_t108_zero = t108 == 0.0
367 if any(mask_t108_zero):
368 t_t108_zero = True
369 t108[mask_t108_zero] = 1.0
370
371 t112 = sqrt(t98 - 2.0 * t99 * t97 + t102 + 2.0 * (t91 * t68 + t95 * t55) * t108)
372
373
374 mask_t112_zero = t112 == 0.0
375 if any(mask_t112_zero):
376 t_t112_zero = True
377 t112[mask_t112_zero] = 1.0
378
379 t97_t99 = t97 + t99
380 t97_nt99 = t97 - t99
381 t113 = t97_nt99 - t112
382 t115 = num_cpmg
383 t116 = power(0.5*(t97_t99 + t112), t115)
384 t118 = 1.0/t112
385 t120 = t97_nt99 + t112
386 t122 = power(0.5*(t97_t99 - t112), t115)
387 t127 = 0.5/t108
388 t120_t122 = t120*t122
389 t139 = 0.5/(k_AB + k_BA) * ((t120_t122 - t113*t116)*t118*k_BA + (t120_t122 - t116*t120)*t127*t113*t118*k_AB)
390
391
392 intensity0 = pA
393 intensity = t139.real * exp(-relax_time * r20)
394
395
396 Mx = intensity / intensity0
397
398
399 back_calc[:] = -inv_relax_time * log(Mx)
400
401
402
403 if t_dw_zero:
404 back_calc[mask_dw_zero.mask] = r20[mask_dw_zero.mask]
405
406
407 if t_t108_zero:
408 back_calc[mask_t108_zero] = 1e100
409
410
411 if t_t112_zero:
412 back_calc[mask_t112_zero] = 1e100
413
414
415
416 if not isfinite(sum(back_calc)):
417
418 fix_invalid(back_calc, copy=False, fill_value=1e100)
419