Source Code for Module lib.dispersion.cr72

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2009 Sebastien Morin                                          # 
  4  # Copyright (C) 2013-2014 Edward d'Auvergne                                   # 
  5  # Copyright (C) 2014 Troels E. Linnet                                         # 
  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  r"""The Carver and Richards (1972) 2-site all time scale exchange U{CR72<http://wiki.nmr-relax.com/CR72>} and U{CR72 full<http://wiki.nmr-relax.com/CR72_full>} models. 
 26   
 27  Description 
 28  =========== 
 29   
 30  This module is for the function, gradient and Hessian of the U{CR72<http://wiki.nmr-relax.com/CR72>} and U{CR72 full<http://wiki.nmr-relax.com/CR72_full>} models. 
 31   
 32   
 33  References 
 34  ========== 
 35   
 36  The model is named after the reference: 
 37   
 38      - Carver, J. P. and Richards, R. E. (1972).  General 2-site solution for chemical exchange produced dependence of T2 upon Carr-Purcell pulse separation.  I{J. Magn. Reson.}, B{6}, 89-105.  (U{DOI: 10.1016/0022-2364(72)90090-X<http://dx.doi.org/10.1016/0022-2364(72)90090-X>}). 
 39   
 40   
 41  Equations 
 42  ========= 
 43   
 44  The equation used is:: 
 45   
 46      R2eff = 1/2 [ R2A0 + R2B0 + kex - 2.nu_cpmg.cosh^-1 (D+.cosh(eta+) - D-.cos(eta-)) ] , 
 47   
 48  where:: 
 49   
 50             1 /        Psi + 2delta_omega^2 \  
 51      D+/- = - | +/-1 + -------------------- | , 
 52             2 \        sqrt(Psi^2 + zeta^2) / 
 53   
 54                             1 
 55      eta+/- = 2^(-3/2) . -------- sqrt(+/-Psi + sqrt(Psi^2 + zeta^2)) , 
 56                          nu_cpmg 
 57   
 58      Psi = (R2A0 - R2B0 - pA.kex + pB.kex)^2 - delta_omega^2 + 4pA.pB.kex^2 , 
 59   
 60      zeta = 2delta_omega (R2A0 - R2B0 - pA.kex + pB.kex). 
 61   
 62  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. 
 63   
 64   
 65  CR72 model 
 66  ---------- 
 67   
 68  Importantly for the implementation of this model, it is assumed that R2A0 and R2B0 are identical.  This simplifies some of the equations to:: 
 69   
 70      R2eff = R20 + kex/2 - nu_cpmg.cosh^-1 (D+.cosh(eta+) - D-.cos(eta-) , 
 71   
 72  where:: 
 73   
 74      Psi = kex^2 - delta_omega^2 , 
 75   
 76      zeta = -2delta_omega (pA.kex - pB.kex). 
 77   
 78   
 79  Links 
 80  ===== 
 81   
 82  More information on the CR72 model can be found in the: 
 83   
 84      - U{relax wiki<http://wiki.nmr-relax.com/CR72>}, 
 85      - U{relax manual<http://www.nmr-relax.com/manual/The_reduced_CR72_2_site_CPMG_model.html>}, 
 86      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#CR72>}. 
 87   
 88  More information on the CR72 full model can be found in the: 
 89   
 90      - U{relax wiki<http://wiki.nmr-relax.com/CR72_full>}, 
 91      - U{relax manual<http://www.nmr-relax.com/manual/The_full_CR72_2_site_CPMG_model.html>}, 
 92      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#CR72_full>}. 
 93  """ 
 94   
 95  # Python module imports. 
 96  from numpy import arccosh, cos, cosh, isfinite, fabs, min, max, multiply, sqrt, subtract, sum 
 97  from numpy.ma import fix_invalid, masked_greater_equal, masked_where 
 98   
 99  # Repetitive calculations (to speed up calculations). 
100  eta_scale = 2.0**(-3.0/2.0) 
101   
102   
103 -def r2eff_CR72(r20a=None, r20a_orig=None, r20b=None, r20b_orig=None, pA=None, dw=None, dw_orig=None, kex=None, cpmg_frqs=None, back_calc=None): 
104      """Calculate the R2eff values for the CR72 model. 
105   
106      See the module docstring for details. 
107   
108   
109      @keyword r20a:          The R20 parameter value of state A (R2 with no exchange). 
110      @type r20a:             numpy float array of rank [NE][NS][NM][NO][ND] 
111      @keyword r20a_orig:     The R20 parameter value of state A (R2 with no exchange). This is only for faster checking of zero value, which result in no exchange. 
112      @type r20a_orig:        numpy float array of rank-1 
113      @keyword r20b:          The R20 parameter value of state B (R2 with no exchange). 
114      @type r20b:             numpy float array of rank [NE][NS][NM][NO][ND] 
115      @keyword r20b_orig:     The R20 parameter value of state B (R2 with no exchange). This is only for faster checking of zero value, which result in no exchange. 
116      @type r20b_orig:        numpy float array of rank-1 
117      @keyword pA:            The population of state A. 
118      @type pA:               float 
119      @keyword dw:            The chemical exchange difference between states A and B in rad/s. 
120      @type dw:               numpy array of rank [NE][NS][NM][NO][ND] 
121      @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. 
122      @type dw_orig:          numpy float array of rank-1 
123      @keyword kex:           The kex parameter value (the exchange rate in rad/s). 
124      @type kex:              float 
125      @keyword cpmg_frqs:     The CPMG nu1 frequencies. 
126      @type cpmg_frqs:        numpy float array of rank [NE][NS][NM][NO][ND] 
127      @keyword back_calc:     The array for holding the back calculated R2eff values.  Each element corresponds to one of the CPMG nu1 frequencies. 
128      @type back_calc:        numpy float array of rank [NE][NS][NM][NO][ND] 
129      """ 
130   
131      # Flag to tell if values should be replaced if max_etapos in cosh function is violated. 
132      t_dw_zero = False 
133      t_max_etapos = False 
134   
135      # Catch parameter values that will result in no exchange, returning flat R2eff = R20 lines (when kex = 0.0, k_AB = 0.0). 
136      # Test if pA or kex is zero. 
137      if kex == 0.0 or pA == 1.0: 
138          back_calc[:] = r20a 
139          return 
140   
141      # 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. 
142      if min(fabs(dw_orig)) == 0.0: 
143          t_dw_zero = True 
144          mask_dw_zero = masked_where(dw == 0.0, dw) 
145   
146      # The B population. 
147      pB = 1.0 - pA 
148   
149      # Repetitive calculations (to speed up calculations). 
150      dw2 = dw**2 
151      r20_kex = (r20a + r20b + kex) / 2.0 
152      k_BA = pA * kex 
153      k_AB = pB * kex 
154   
155      # The Psi and zeta values. 
156      if sum(r20a_orig - r20b_orig) != 0.0: 
157          fact = r20a - r20b - k_BA + k_AB 
158          Psi = fact**2 - dw2 + 4.0*k_BA*k_AB 
159          zeta = 2.0*dw * fact 
160      else: 
161          Psi = kex**2 - dw2 
162          zeta = -2.0*dw * (k_BA - k_AB) 
163   
164      # More repetitive calculations. 
165      sqrt_psi2_zeta2 = sqrt(Psi**2 + zeta**2) 
166   
167      # The D+/- values. 
168      D_part = (0.5*Psi + dw2) / sqrt_psi2_zeta2 
169      Dpos = 0.5 + D_part 
170      Dneg = -0.5 + D_part 
171   
172      # Partial eta+/- values. 
173      eta_fact = eta_scale / cpmg_frqs 
174      etapos = eta_fact * sqrt(Psi + sqrt_psi2_zeta2) 
175      etaneg = eta_fact * sqrt(-Psi + sqrt_psi2_zeta2) 
176   
177      # Catch math domain error of cosh(val > 710). 
178      # This is when etapos > 710. 
179      if max(etapos) > 700: 
180          t_max_etapos = True 
181          mask_max_etapos = masked_greater_equal(etapos, 700.0) 
182          # To prevent math errors, set etapos to 1. 
183          etapos[mask_max_etapos.mask] = 1.0 
184   
185      # The arccosh argument - catch invalid values. 
186      fact = Dpos * cosh(etapos) - Dneg * cos(etaneg) 
187      if min(fact) < 1.0: 
188          back_calc[:] = r20_kex 
189          return 
190   
191      # Calculate R2eff. This uses the temporary buffer and fill directly to back_calc. 
192      multiply(cpmg_frqs, arccosh(fact), out=back_calc) 
193      subtract(r20_kex, back_calc, out=back_calc) 
194   
195      # Replace data in array. 
196      # If dw is zero. 
197      if t_dw_zero: 
198          back_calc[mask_dw_zero.mask] = r20a[mask_dw_zero.mask] 
199   
200      # If eta_pos above 700. 
201      if t_max_etapos: 
202          back_calc[mask_max_etapos.mask] = r20a[mask_max_etapos.mask] 
203   
204      # Catch errors, taking a sum over array is the fastest way to check for 
205      # +/- inf (infinity) and nan (not a number). 
206      if not isfinite(sum(back_calc)): 
207          # Replaces nan, inf, etc. with fill value. 
208          fix_invalid(back_calc, copy=False, fill_value=1e100) 
209