Source Code for Module lib.dispersion.b14

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