Source Code for Module lib.dispersion.tsmfk01

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2013-2014 Edward d'Auvergne                                   # 
  4  # Copyright (C) 2013-2014 Troels E. Linnet                                    # 
  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 Tollinger et al (2001) 2-site very-slow exchange U{TSMFK01<http://wiki.nmr-relax.com/TSMFK01>} model. 
 25   
 26  Description 
 27  =========== 
 28   
 29  Applicable in the limit of slow exchange, range of microsecond to second time scale, when |R2A-R2B| << k_AB, kB << 1/tau_CP.  R20A is the transverse relaxation rate of site A in the absence of exchange.  2*tau_CP is is the time between successive 180 degree pulses. 
 30   
 31  This module is for the function, gradient and Hessian of the U{TSMFK01<http://wiki.nmr-relax.com/TSMFK01>} model 
 32   
 33   
 34  References 
 35  ========== 
 36   
 37  The model is named after the reference: 
 38   
 39      - Tollinger, M., Skrynnikov, N. R., Mulder, F. A. A., Forman-Kay, J. D. and Kay, L. E. (2001).  Slow Dynamics in Folded and Unfolded States of an SH3 Domain, I{J. Am. Chem. Soc.}, B{123} (46) (U{DOI: 10.1021/ja011300z<http://dx.doi.org/10.1021/ja011300z>}). 
 40   
 41   
 42  Equations 
 43  ========= 
 44   
 45  The equation used is:: 
 46   
 47                                     sin(delta_omega * tau_CP) 
 48      R2Aeff = R20A + k_AB - k_AB * -------------------------  , 
 49                                     delta_omega * tau_CP 
 50   
 51  where:: 
 52   
 53      tau_CP = 1.0/(4*nu_cpmg) , 
 54   
 55  R20A is the transverse relaxation rate of site A in the absence of exchange, 2*tau_CP is is the time between successive 180 deg. pulses, k_AB is the forward chemical exchange rate constant, delta_omega is the chemical shift difference between the two states. 
 56   
 57   
 58  Links 
 59  ===== 
 60   
 61  More information on the TSMFK01 model can be found in the: 
 62   
 63      - U{relax wiki<http://wiki.nmr-relax.com/TSMFK01>}, 
 64      - U{relax manual<http://www.nmr-relax.com/manual/The_TSMFK01_2_site_CPMG_model.html>}, 
 65      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#TSMFK01>}. 
 66  """ 
 67   
 68  # Python module imports. 
 69  from numpy import fabs, min, sin, isfinite, sum 
 70  from numpy.ma import fix_invalid, masked_where 
 71   
 72   
 73 -def r2eff_TSMFK01(r20a=None, dw=None, dw_orig=None, k_AB=None, tcp=None, back_calc=None): 
 74      """Calculate the R2eff values for the TSMFK01 model. 
 75   
 76      See the module docstring for details. 
 77   
 78   
 79      @keyword r20a:          The R20 parameter value of state A (R2 with no exchange). 
 80      @type r20a:             numpy float array of rank [NE][NS][NM][NO][ND] 
 81      @keyword dw:            The chemical exchange difference between states A and B in rad/s. 
 82      @type dw:               numpy float array of rank [NE][NS][NM][NO][ND] 
 83      @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. 
 84      @type dw_orig:          numpy float array of rank-1 
 85      @keyword k_AB:          The k_AB parameter value (the forward exchange rate in rad/s). 
 86      @type k_AB:             float 
 87      @keyword tcp:           The tau_CPMG times (1 / 4.nu1). 
 88      @type tcp:              numpy float array of rank [NE][NS][NM][NO][ND] 
 89      @keyword back_calc:     The array for holding the back calculated R2eff values.  Each element corresponds to one of the CPMG nu1 frequencies. 
 90      @type back_calc:        numpy float array of rank [NE][NS][NM][NO][ND] 
 91      """ 
 92   
 93      # Flag to tell if values should be replaced if max_etapos in cosh function is violated. 
 94      t_dw_zero = False 
 95   
 96      # Catch parameter values that will result in no exchange, returning flat R2eff = R20 lines (when kex = 0.0, k_AB = 0.0). 
 97      # Test if k_AB is zero. 
 98      if k_AB == 0.0: 
 99          back_calc[:] = r20a 
100          return 
101   
102      # 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. 
103      if min(fabs(dw_orig)) == 0.0: 
104          t_dw_zero = True 
105          mask_dw_zero = masked_where(dw == 0.0, dw) 
106   
107      # Denominator. 
108      denom = dw * tcp 
109   
110      # The numerator. 
111      numer = sin(denom) 
112   
113      # Catch zeros (to avoid pointless mathematical operations). 
114      # This will result in no exchange, returning flat lines. 
115      if min(fabs(numer)) == 0.0: 
116          # Calculate R2eff for forward. 
117          back_calc[:] = r20a + k_AB 
118      else: 
119          # Calculate R2eff. 
120          back_calc[:] = r20a + k_AB - k_AB * numer / denom 
121   
122      # Replace data in array. 
123      # If dw is zero. 
124      if t_dw_zero: 
125          back_calc[mask_dw_zero.mask] = r20a[mask_dw_zero.mask] 
126   
127      # Catch errors, taking a sum over array is the fastest way to check for 
128      # +/- inf (infinity) and nan (not a number). 
129      if not isfinite(sum(back_calc)): 
130          # Replaces nan, inf, etc. with fill value. 
131          fix_invalid(back_calc, copy=False, fill_value=1e100) 
132