Source Code for Module lib.dispersion.it99

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2009 Sebastien Morin                                          # 
  4  # Copyright (C) 2013-2014 Edward d'Auvergne                                   # 
  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 Ishima and Torchia (1999) 2-site all time scale exchange (with pA >> pB) U{IT99<http://wiki.nmr-relax.com/IT99>} model. 
 25   
 26  Description 
 27  =========== 
 28   
 29  This module is for the function, gradient and Hessian of the U{IT99<http://wiki.nmr-relax.com/IT99>} model. 
 30   
 31   
 32  References 
 33  ========== 
 34   
 35  The model is named after the reference: 
 36   
 37      - 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>}). 
 38   
 39   
 40  Equations 
 41  ========= 
 42   
 43  The equation used is:: 
 44   
 45                phi_ex * tex 
 46      Rex ~= ------------------- , 
 47             1 + omega_a^2*tex^2 
 48   
 49      phi_ex = pA * pB * delta_omega^2 , 
 50   
 51      omega_a^2 = sqrt(omega_1eff^4 + pA^2*delta_omega^4) , 
 52   
 53      R2eff = R20 + Rex , 
 54   
 55  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:: 
 56   
 57      omega_1eff = 4*sqrt(3) * nu_cpmg 
 58   
 59  and therefore:: 
 60   
 61      omega_1eff^4 = 2304 * nu_cpmg^4 
 62   
 63   
 64  Links 
 65  ===== 
 66   
 67  More information on the IT99 model can be found in the: 
 68   
 69      - U{relax wiki<http://wiki.nmr-relax.com/IT99>}, 
 70      - U{relax manual<http://www.nmr-relax.com/manual/IT99_2_site_CPMG_model.html>}, 
 71      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#IT99>}. 
 72   
 73  """ 
 74   
 75  # Python module imports. 
 76  from numpy import array, isfinite, sqrt, sum 
 77   
 78   
 79 -def r2eff_IT99(r20=None, pA=None, pB=None, dw=None, tex=None, cpmg_frqs=None, back_calc=None, num_points=None): 
 80      """Calculate the R2eff values for the IT99 model. 
 81   
 82      See the module docstring for details. 
 83   
 84   
 85      @keyword r20:           The R20 parameter value (R2 with no exchange). 
 86      @type r20:              float 
 87      @keyword pA:            The population of state A. 
 88      @type pA:               float 
 89      @keyword pB:            The population of state B. 
 90      @type pB:               float 
 91      @keyword dw:            The chemical exchange difference between states A and B in rad/s. 
 92      @type dw:               float 
 93      @keyword tex:           The tex parameter value (the time of exchange in s/rad). 
 94      @type tex:              float 
 95      @keyword cpmg_frqs:     The CPMG nu1 frequencies. 
 96      @type cpmg_frqs:        numpy rank-1 float array 
 97      @type tcp:              numpy rank-1 float array 
 98      @keyword back_calc:     The array for holding the back calculated R2eff values.  Each element corresponds to one of the CPMG nu1 frequencies. 
 99      @type back_calc:        numpy rank-1 float array 
100      @keyword num_points:    The number of points on the dispersion curve, equal to the length of the cpmg_frqs and back_calc arguments. 
101      @type num_points:       int 
102      """ 
103   
104      # Repetitive calculations (to speed up calculations). 
105      dw2 = dw**2 
106      tex2 = tex**2 
107      padw2 = pA * dw2 
108      pa2dw4 = padw2**2 
109   
110      # The numerator. 
111      numer = padw2 * pB * tex 
112   
113      # Catch zeros (to avoid pointless mathematical operations). 
114      # This will result in no exchange, returning flat lines. 
115      if numer == 0.0: 
116          back_calc[:] = array([r20]*num_points) 
117          return 
118   
119      # The effective rotating frame field. 
120      omega_1eff4 = 2304.0 * cpmg_frqs**4 
121   
122      # Denominator. 
123      omega_a2 = sqrt(omega_1eff4 + pa2dw4) 
124      denom = 1.0 + omega_a2 * tex2 
125   
126      # R2eff calculation. 
127      R2eff = r20 + numer / denom 
128   
129      # Catch errors, taking a sum over array is the fastest way to check for 
130      # +/- inf (infinity) and nan (not a number). 
131      if not isfinite(sum(R2eff)): 
132          R2eff = array([1e100]*num_points) 
133   
134      back_calc[:] = R2eff 
135