Source Code for Module lib.dispersion.it99

  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  """The Ishima and Torchia (1999) 2-site all time scale exchange (with pA >> pB) U{IT99<http://wiki.nmr-relax.com/IT99>} model. 
 26   
 27  Description 
 28  =========== 
 29   
 30  This module is for the function, gradient and Hessian of the U{IT99<http://wiki.nmr-relax.com/IT99>} model. 
 31   
 32   
 33  References 
 34  ========== 
 35   
 36  The model is named after the reference: 
 37   
 38      - Ishima R. and Torchia D.A. (1999).  Estimating the time scale of chemical exchange of proteins from measurements of transverse relaxation rates in solution.  I{J. Biomol. NMR}, B{14}, 369-372.  (U{DOI: 10.1023/A:1008324025406<http://dx.doi.org/10.1023/A:1008324025406>}). 
 39   
 40   
 41  Equations 
 42  ========= 
 43   
 44  The equation used is:: 
 45   
 46                phi_ex * tex 
 47      Rex ~= ------------------- , 
 48             1 + omega_a^2*tex^2 
 49   
 50      phi_ex = pA * pB * delta_omega^2 , 
 51   
 52      omega_a^2 = sqrt(omega_1eff^4 + pA^2*delta_omega^4) , 
 53   
 54      R2eff = R20 + Rex , 
 55   
 56  where tex = 1/(2kex), 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.  The effective rotating frame field for a CPMG-type experiment is given by:: 
 57   
 58      omega_1eff = 4*sqrt(3) * nu_cpmg 
 59   
 60  and therefore:: 
 61   
 62      omega_1eff^4 = 2304 * nu_cpmg^4 
 63   
 64   
 65  Links 
 66  ===== 
 67   
 68  More information on the IT99 model can be found in the: 
 69   
 70      - U{relax wiki<http://wiki.nmr-relax.com/IT99>}, 
 71      - U{relax manual<http://www.nmr-relax.com/manual/The_IT99_2_site_CPMG_model.html>}, 
 72      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#IT99>}. 
 73   
 74  """ 
 75   
 76  # Python module imports. 
 77  from numpy import isfinite, fabs, sqrt, sum 
 78  from numpy.ma import fix_invalid, masked_where 
 79   
 80   
 81 -def r2eff_IT99(r20=None, pA=None, dw=None, dw_orig=None, tex=None, cpmg_frqs=None, back_calc=None): 
 82      """Calculate the R2eff values for the IT99 model. 
 83   
 84      See the module docstring for details. 
 85   
 86   
 87      @keyword r20:           The R20 parameter value (R2 with no exchange). 
 88      @type r20:              numpy float array of rank [NE][NS][NM][NO][ND] 
 89      @keyword pA:            The population of state A. 
 90      @type pA:               float 
 91      @keyword dw:            The chemical exchange difference between states A and B in rad/s. 
 92      @type dw:               numpy float array of rank [NE][NS][NM][NO][ND] 
 93      @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. 
 94      @type dw_orig:          numpy float array of rank-1 
 95      @keyword tex:           The tex parameter value (the time of exchange in s/rad). 
 96      @type tex:              float 
 97      @keyword cpmg_frqs:     The CPMG nu1 frequencies. 
 98      @type cpmg_frqs:        numpy float array of rank [NE][NS][NM][NO][ND] 
 99      @keyword back_calc:     The array for holding the back calculated R2eff values.  Each element corresponds to one of the CPMG nu1 frequencies. 
100      @type back_calc:        numpy float array of rank [NE][NS][NM][NO][ND] 
101      """ 
102   
103      # Flag to tell if values should be replaced if numer is zero. 
104      t_dw_zero = False 
105   
106      # Catch divide with zeros (to avoid pointless mathematical operations). 
107      if tex == 0.0 or pA == 1.0: 
108          back_calc[:] = r20 
109          return 
110   
111      # 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. 
112      if min(fabs(dw_orig)) == 0.0: 
113          t_dw_zero = True 
114          mask_dw_zero = masked_where(dw == 0.0, dw) 
115   
116      # Parameter conversions. 
117      pB = 1.0 - pA 
118   
119      # Repetitive calculations (to speed up calculations). 
120      padw2 = pA * dw**2 
121      pa2dw4 = padw2**2 
122   
123      # The numerator. 
124      numer = padw2 * pB * tex 
125   
126      # The effective rotating frame field. 
127      omega_1eff4 = 2304.0 * cpmg_frqs**4 
128   
129      # Denominator. 
130      omega_a2 = sqrt(omega_1eff4 + pa2dw4) 
131      denom = 1.0 + omega_a2 * tex**2 
132   
133      # R2eff calculation. 
134      back_calc[:] = r20 + numer / denom 
135   
136      # Replace data in array. 
137      # If dw is zero. 
138      if t_dw_zero: 
139          back_calc[mask_dw_zero.mask] = r20[mask_dw_zero.mask] 
140   
141      # Catch errors, taking a sum over array is the fastest way to check for 
142      # +/- inf (infinity) and nan (not a number). 
143      if not isfinite(sum(back_calc)): 
144          # Replaces nan, inf, etc. with fill value. 
145          fix_invalid(back_calc, copy=False, fill_value=1e100) 
146