Package lib :: Package dispersion :: Module ns_cpmg_2site_expanded
[hide private]
[frames] | no frames]

Source Code for Module lib.dispersion.ns_cpmg_2site_expanded

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2000-2001 Nikolai Skrynnikov                                  # 
  4  # Copyright (C) 2000-2001 Martin Tollinger                                    # 
  5  # Copyright (C) 2010-2013 Paul Schanda                                        # 
  6  # Copyright (C) 2013 Mathilde Lescanne                                        # 
  7  # Copyright (C) 2013 Dominique Marion                                         # 
  8  # Copyright (C) 2013-2014 Edward d'Auvergne                                   # 
  9  # Copyright (C) 2014 Troels E. Linnet                                         # 
 10  #                                                                             # 
 11  # This file is part of the program relax (http://www.nmr-relax.com).          # 
 12  #                                                                             # 
 13  # This program is free software: you can redistribute it and/or modify        # 
 14  # it under the terms of the GNU General Public License as published by        # 
 15  # the Free Software Foundation, either version 3 of the License, or           # 
 16  # (at your option) any later version.                                         # 
 17  #                                                                             # 
 18  # This program is distributed in the hope that it will be useful,             # 
 19  # but WITHOUT ANY WARRANTY without even the implied warranty of               # 
 20  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 21  # GNU General Public License for more details.                                # 
 22  #                                                                             # 
 23  # You should have received a copy of the GNU General Public License           # 
 24  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 25  #                                                                             # 
 26  ############################################################################### 
 27   
 28  # Module docstring. 
 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  # Python module imports. 
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 # Flag to tell if values should be replaced if math function is violated. 283 t_dw_zero = False 284 t_t108_zero = False 285 t_t112_zero = False 286 287 # Catch parameter values that will result in no exchange, returning flat R2eff = R20 lines (when kex = 0.0, k_AB = 0.0). 288 if pA == 1.0 or kex == 0.0: 289 back_calc[:] = r20 290 return 291 292 # Test if dw is zero. Create a mask for the affected spins to replace these with R20 at the end of the calculationWait for replacement, since this is spin specific. 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 # Once off parameter conversions. 298 pB = 1.0 - pA 299 k_BA = pA * kex 300 k_AB = pB * kex 301 302 # Repetitive calculations. 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 # The expansion factors (in numpy array form for all dispersion points). 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 # If t108 is zero. 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 # If t112 is zero. 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 # Calculate the initial and final peak intensities. 392 intensity0 = pA 393 intensity = t139.real * exp(-relax_time * r20) 394 395 # The magnetisation vector. 396 Mx = intensity / intensity0 397 398 # Calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential, and store it for each dispersion point. 399 back_calc[:] = -inv_relax_time * log(Mx) 400 401 # Replace data in array. 402 # If dw is zero. 403 if t_dw_zero: 404 back_calc[mask_dw_zero.mask] = r20[mask_dw_zero.mask] 405 406 # If t108 is zero. 407 if t_t108_zero: 408 back_calc[mask_t108_zero] = 1e100 409 410 # If t112 is zero. 411 if t_t112_zero: 412 back_calc[mask_t112_zero] = 1e100 413 414 # Catch errors, taking a sum over array is the fastest way to check for 415 # +/- inf (infinity) and nan (not a number). 416 if not isfinite(sum(back_calc)): 417 # Replaces nan, inf, etc. with fill value. 418 fix_invalid(back_calc, copy=False, fill_value=1e100)
419