Source Code for Module lib.dispersion.b14

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2014 Troels E. Linnet                                         # 
  4  # Copyright (C) 2014 Andrew Baldwin                                           # 
  5  #                                                                             # 
  6  # This file is part of the program relax (http://www.nmr-relax.com).          # 
  7  #                                                                             # 
  8  # This program is free software: you can redistribute it and/or modify        # 
  9  # it under the terms of the GNU General Public License as published by        # 
 10  # the Free Software Foundation, either version 3 of the License, or           # 
 11  # (at your option) any later version.                                         # 
 12  #                                                                             # 
 13  # This program is distributed in the hope that it will be useful,             # 
 14  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 15  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 16  # GNU General Public License for more details.                                # 
 17  #                                                                             # 
 18  # You should have received a copy of the GNU General Public License           # 
 19  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 20  #                                                                             # 
 21  ############################################################################### 
 22   
 23  # Module docstring. 
 24  """The Baldwin (2014) 2-site exact solution model for all time scales U{B14<http://wiki.nmr-relax.com/B14>}. 
 25   
 26  Description 
 27  =========== 
 28   
 29  This module is for the function, gradient and Hessian of the U{B14<http://wiki.nmr-relax.com/B14>} model. 
 30   
 31   
 32  References 
 33  ========== 
 34   
 35  The model is named after the reference: 
 36   
 37      - Andrew J. Baldwin (2014).  An exact solution for R2,eff in CPMG experiments in the case of two site chemical exchange.  I{J. Magn. Reson.}.  (U{DOI: 10.1016/j.jmr.2014.02.023 <http://dx.doi.org/10.1016/j.jmr.2014.02.023>}). 
 38   
 39   
 40  Equations 
 41  ========= 
 42   
 43  The equation used is:: 
 44   
 45              R2A0 + R2B0 + kex      Ncyc                      1      ( 1+y            1-y                          ) 
 46      R2eff = ------------------ -  ------ * cosh^-1 * v1c - ------ ln( --- + ------------------ * (v2 + 2*kAB*pD ) ) , 
 47                    2                Trel                     Trel    (  2    2*sqrt(v1c^2 -1 )                     ) 
 48   
 49                              1      ( 1+y            1-y                          ) 
 50            = R2eff(CR72) - ------ ln( --- + ------------------ * (v2 + 2*kAB*pD ) ) , 
 51                             Trel    (  2    2*sqrt(v1c^2 -1 )                     ) 
 52   
 53  Which have these following definitions:: 
 54   
 55      v1c = F0 * cosh(tauCP * E0)- F2 * cosh(tauCP * E2) , 
 56      v1s = F0 * sinh(tauCP * E0)- F2 * sinh(tauCP * E2) , 
 57      v2*N = v1s * (OB-OA) + 4OB * F1^a * sinh(tauCP * E1) , 
 58      pD N = v1s + (F1^a + F1^b) * sinh(tauCP * E1) , 
 59      v3 = ( v2^2 + 4 * kBA * kAB * pD^2 )^1/2 , 
 60      y = ( (v1c-v3)/(v1c+v3) )^NCYC , 
 61   
 62  Note, E2 is complex. If |x| denotes the complex modulus:: 
 63   
 64      cosh(tauCP * E2) = cos(tauCP * |E2|) , 
 65      sinh(tauCP * E2) = i sin(tauCP * |E2|) , 
 66   
 67  The term pD is based on product of the off diagonal elements in the CPMG propagator (Supplementary Section 3). 
 68   
 69  It is interesting to consider the region of validity of the Carver Richards result.  The two results are equal when the correction is zero, which is true when:: 
 70   
 71      sqrt(v1c^2-1) ~ v2 + 2*kAB * pD , 
 72   
 73  This occurs when 2*kAB * pD tends to zero, and so v2=v3.  Setting "kAB * pD" to zero, amounts to neglecting magnetisation that starts on the ground state ensemble and end on the excited state ensemble and vice versa.  This will be a good approximation when pA >> p_B. 
 74   
 75  In practise, significant deviations from the Carver Richards equation can be incurred if pB > 1 %.  Incorporation of the correction term into equation (50), results in an improved description of the CPMG experiment over the Carver Richards equation. 
 76   
 77  kex is the chemical exchange rate constant, pA and pB are the populations of states A and B, and delta_omega is the chemical shift difference between the two states in ppm. 
 78   
 79   
 80  Links 
 81  ===== 
 82   
 83  More information on the B14 model can be found in the: 
 84   
 85      - U{relax wiki<http://wiki.nmr-relax.com/B14>}, 
 86      - U{relax manual<http://www.nmr-relax.com/manual/The_reduced_B14_2_site_CPMG_model.html>}, 
 87      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#B14>}. 
 88   
 89  More information on the B14 full model can be found in the: 
 90   
 91      - U{relax wiki<http://wiki.nmr-relax.com/B14_full>}, 
 92      - U{relax manual<http://www.nmr-relax.com/manual/The_full_B14_2_site_CPMG_model.html>}, 
 93      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#B14_full>}. 
 94   
 95   
 96  Comparison to CR72 model 
 97  ======================== 
 98   
 99  Comparison to CR72 model can be found in the: 
100   
101      - U{relax wiki<http://wiki.nmr-relax.com/CR72>}, 
102      - U{relax manual<http://www.nmr-relax.com/manual/The_reduced_CR72_2_site_CPMG_model.html>}, 
103      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#CR72>}. 
104   
105  Comparison to CR72 full model can be found in the: 
106   
107      - U{relax wiki<http://wiki.nmr-relax.com/CR72_full>}, 
108      - U{relax manual<http://www.nmr-relax.com/manual/The_full_CR72_2_site_CPMG_model.html>}, 
109      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#CR72_full>}. 
110  """ 
111   
112  # Python module imports. 
113  from numpy import any, arccosh, arctan2, cos, cosh, fabs, isfinite, log, max, min, power, sin, sinh, sqrt, sum 
114  from numpy.ma import fix_invalid, masked_greater_equal, masked_where 
115   
116  # Repetitive calculations (to speed up calculations). 
117  g_fact = 1/sqrt(2) 
118   
119   
120 -def r2eff_B14(r20a=None, r20b=None, pA=None, dw=None, dw_orig=None, kex=None, ncyc=None, inv_tcpmg=None, tcp=None, back_calc=None): 
121      """Calculate the R2eff values for the CR72 model. 
122   
123      See the module docstring for details. 
124   
125   
126      @keyword r20a:          The R20 parameter value of state A (R2 with no exchange). 
127      @type r20a:             numpy float array of rank [NE][NS][NM][NO][ND] 
128      @keyword r20b:          The R20 parameter value of state B (R2 with no exchange). 
129      @type r20b:             numpy float array of rank [NE][NS][NM][NO][ND] 
130      @keyword pA:            The population of state A. 
131      @type pA:               float 
132      @keyword dw:            The chemical exchange difference between states A and B in rad/s. 
133      @type dw:               numpy float array of rank [NE][NS][NM][NO][ND] 
134      @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. 
135      @type dw_orig:          numpy float array of rank-1 
136      @keyword kex:           The kex parameter value (the exchange rate in rad/s). 
137      @type kex:              float 
138      @keyword ncyc:          The matrix exponential power array. The number of CPMG blocks. 
139      @type ncyc:             numpy int16 array of rank [NE][NS][NM][NO][ND] 
140      @keyword inv_tcpmg:     The inverse of the total duration of the CPMG element (in inverse seconds). 
141      @type inv_tcpmg:        numpy float array of rank [NE][NS][NM][NO][ND] 
142      @keyword tcp:           The tau_CPMG times (1 / 4.nu1). 
143      @type tcp:              numpy float array of rank [NE][NS][NM][NO][ND] 
144      @keyword back_calc:     The array for holding the back calculated R2eff values.  Each element corresponds to one of the CPMG nu1 frequencies. 
145      @type back_calc:        numpy float array of rank [NE][NS][NM][NO][ND] 
146      """ 
147   
148      # Flag to tell if values should be replaced if math function is violated. 
149      t_dw_zero = False 
150      t_max_e = False 
151      t_v3_N_zero = False 
152      t_log_tog_neg = False 
153      t_v1c_less_one = False 
154   
155      # Catch parameter values that will result in no exchange, returning flat R2eff = R20 lines (when kex = 0.0, k_AB = 0.0). 
156      # Test if pA or kex is zero. 
157      if kex == 0.0 or pA == 1.0: 
158          back_calc[:] = r20a 
159          return 
160   
161      # 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. 
162      if min(fabs(dw_orig)) == 0.0: 
163          t_dw_zero = True 
164          mask_dw_zero = masked_where(dw == 0.0, dw) 
165   
166      # Parameter conversions. 
167      pB = 1.0 - pA 
168      k_BA = pA * kex 
169      k_AB = pB * kex 
170   
171      # Repetitive calculations (to speed up calculations). 
172      deltaR2 = r20a - r20b 
173      dw2 = dw**2 
174      two_tcp = 2.0 * tcp 
175   
176      # The Carver and Richards (1972) alpha_minus short notation. 
177      alpha_m = deltaR2 + k_AB - k_BA 
178      zeta = 2.0 * dw * alpha_m 
179      Psi = alpha_m**2 + 4.0 * k_BA * k_AB - dw2 
180   
181      # Get the real and imaginary components of the exchange induced shift. 
182      # Trigonometric functions faster than square roots. 
183      quad_zeta2_Psi2 = (zeta**2 + Psi**2)**0.25 
184      fact = 0.5 * arctan2(-zeta, Psi) 
185      g3 = cos(fact) * quad_zeta2_Psi2 
186      g4 = sin(fact) * quad_zeta2_Psi2 
187   
188      # Repetitive calculations (to speed up calculations). 
189      g32 = g3**2 
190      g42 = g4**2 
191   
192      # Time independent factors. 
193      # N = oG + oE. 
194      N = g3 + g4*1j 
195   
196      NNc = g32 + g42 
197   
198      # F0. 
199      F0 = (dw2 + g32) / NNc 
200   
201      # F2. 
202      F2 = (dw2 - g42) / NNc 
203   
204      # t1 = (-dw + g4) * (complex(-dw, -g3)) / NNc #t1. 
205   
206      # t2. 
207      F1b = (dw + g4) * (dw - g3*1j) / NNc 
208   
209      # t1 + t2. 
210      F1a_plus_b = (2. * dw2 + zeta*1j) / NNc 
211   
212      # Derived from relaxation. 
213      # E0 = -2.0 * tcp * (F00R - f11R). 
214      E0 = two_tcp * g3 
215   
216      # Catch math domain error of sinh(val > 710). 
217      # This is when E0 > 710. 
218      if max(E0) > 700: 
219          t_max_e = True 
220          mask_max_e = masked_greater_equal(E0, 700.0) 
221          # To prevent math errors, set e_zero to 1. 
222          E0[mask_max_e.mask] = 1.0 
223   
224      # Derived from chemical shifts  #E2 = complex(0, -2.0 * tcp * (F00I - f11I)). 
225      E2 = two_tcp * g4 
226   
227      # Mixed term (complex) (E0 - iE2)/2. 
228      E1 = (g3 - g4*1j) * tcp 
229   
230      # Complex. 
231      v1s = F0 * sinh(E0) - F2 * sin(E2)*1j 
232   
233      # -2 * oG * t2. 
234      v4 = F1b * (-alpha_m - g3 ) + F1b * (dw - g4)*1j 
235   
236      # Complex. 
237      ex1c = sinh(E1) 
238   
239      # Off diagonal common factor. sinh fuctions. 
240      v5 = (-deltaR2 + kex + dw*1j) * v1s - 2. * (v4 + k_AB * F1a_plus_b) * ex1c 
241   
242      # Real. The v_1c in paper. 
243      v1c = F0 * cosh(E0) - F2 * cos(E2) 
244   
245      # Catch math domain error of sqrt of negative. 
246      # This is when v1c is less than 1. 
247      mask_v1c_less_one = v1c < 1.0 
248      if any(mask_v1c_less_one): 
249          t_v1c_less_one = True 
250          v1c[mask_v1c_less_one] = 1.0 
251   
252      # Exact result for v2v3. 
253      v3 = sqrt(v1c**2 - 1.) 
254   
255      y = power( (v1c - v3) / (v1c + v3), ncyc) 
256   
257      Tog_div = 2. * v3 * N 
258   
259      # Catch math domain error of division with 0. 
260      # This is when Tog_div is zero. 
261      mask_v3_N_zero = Tog_div == 0.0 
262      if any(mask_v3_N_zero): 
263          t_v3_N_zero = True 
264          Tog_div[mask_v3_N_zero] = 1.0 
265   
266      Tog = 0.5 * (1. + y) + (1. - y) * v5 / Tog_div 
267   
268      ## -1/Trel * log(LpreDyn). 
269      # Rpre = (r20a + r20b + kex) / 2.0 
270   
271      ## Carver and Richards (1972) 
272      # R2eff_CR72 = Rpre - inv_tcpmg * ncyc *  arccosh(v1c.real) 
273   
274      ## Baldwin final. 
275      # Estimate R2eff. relax_time = Trel = 1/inv_tcpmg. 
276      # R2eff = R2eff_CR72 - inv_tcpmg * log(Tog.real) 
277   
278      # Catch math domain error of log of negative. 
279      # This is when Tog.real is negative. 
280      mask_log_tog_neg = Tog.real < 0.0 
281      if any(mask_log_tog_neg): 
282          t_log_tog_neg = True 
283          Tog.real[mask_log_tog_neg] = 1.0 
284   
285      # Fastest calculation. 
286      back_calc[:] = (r20a + r20b + kex) / 2.0 - inv_tcpmg * ( ncyc * arccosh(v1c.real) + log(Tog.real) ) 
287   
288      # Replace data in array. 
289      # If dw is zero. 
290      if t_dw_zero: 
291          back_calc[mask_dw_zero.mask] = r20a[mask_dw_zero.mask] 
292   
293      # If E0 is above 700. 
294      if t_max_e: 
295          back_calc[mask_max_e.mask] = r20a[mask_max_e.mask] 
296   
297      # If v1c is less than 1. 
298      if t_v1c_less_one: 
299          back_calc[mask_v1c_less_one] = 1e100 
300   
301      # If Tog_div is zero. 
302      if t_v3_N_zero: 
303          back_calc[mask_v3_N_zero] = 1e100 
304   
305      # If Tog.real is negative. 
306      if t_log_tog_neg: 
307          back_calc[mask_log_tog_neg] = 1e100 
308   
309      # Catch errors, taking a sum over array is the fastest way to check for 
310      # +/- inf (infinity) and nan (not a number). 
311      if not isfinite(sum(back_calc)): 
312          # Replaces nan, inf, etc. with fill value. 
313          fix_invalid(back_calc, copy=False, fill_value=1e100) 
314