Source Code for Module lib.dispersion.dpl94

  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 Davis, Perlman and London (1994) 2-site fast exchange R1rho U{DPL94<http://wiki.nmr-relax.com/DPL94>} model. 
 26   
 27  Description 
 28  =========== 
 29   
 30  This module is for the function, gradient and Hessian of the U{DPL94<http://wiki.nmr-relax.com/DPL94>} model. 
 31   
 32   
 33  References 
 34  ========== 
 35   
 36  The model is named after the reference: 
 37   
 38      - Davis, D. G., Perlman, M. E. and London, R. E. (1994).  Direct measurements of the dissociation-rate constant for inhibitor-enzyme complexes via the T1rho and T2 (CPMG) methods.  I{J. Magn. Reson.}, Series B, B{104}, 266-275.  (U{DOI: 10.1006/jmrb.1994.1084<http://dx.doi.org/10.1006/jmrb.1994.1084>}) 
 39   
 40  Equations 
 41  ========= 
 42   
 43  The equation used is:: 
 44   
 45                                                                        phi_ex * kex 
 46      R1rho = R1.cos^2(theta) + R1rho'.sin^2(theta) + sin^2(theta) * ------------------ , 
 47                                                                     kex^2 + omega_sl^2 
 48   
 49  where theta is the rotating frame tilt angle, and:: 
 50   
 51      phi_ex = pA * pB * delta_omega^2 , 
 52   
 53  kex is the chemical exchange rate constant, pA and pB are the populations of states A and B, delta_omega is the chemical shift difference between the two states, and omega_sl is the spin-lock field strength. 
 54   
 55   
 56  Links 
 57  ===== 
 58   
 59  More information on the DPL94 model can be found in the: 
 60   
 61      - U{relax wiki<http://wiki.nmr-relax.com/DPL94>}, 
 62      - U{relax manual<http://www.nmr-relax.com/manual/The_DPL94_2_site_fast_exchange_R1_rho_model.html>}, 
 63      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#DPL94>}. 
 64  """ 
 65   
 66  # Python module imports. 
 67  from numpy import any, cos, isfinite, min, sin, sum 
 68  from numpy.ma import fix_invalid, masked_where 
 69   
 70   
 71 -def r1rho_DPL94(r1rho_prime=None, phi_ex=None, kex=None, theta=None, R1=0.0, spin_lock_fields2=None, back_calc=None): 
 72      """Calculate the R1rho values for the DPL94 model. 
 73   
 74      See the module docstring for details. 
 75   
 76   
 77      @keyword r1rho_prime:       The R1rho_prime parameter value (R1rho with no exchange). 
 78      @type r1rho_prime:          numpy float array of rank [NE][NS][NM][NO][ND] 
 79      @keyword phi_ex:            The phi_ex parameter value (pA * pB * delta_omega^2). 
 80      @type phi_ex:               numpy float array of rank [NE][NS][NM][NO][ND] 
 81      @keyword kex:               The kex parameter value (the exchange rate in rad/s). 
 82      @type kex:                  float 
 83      @keyword theta:             The rotating frame tilt angles for each dispersion point. 
 84      @type theta:                numpy float array of rank [NE][NS][NM][NO][ND] 
 85      @keyword R1:                The R1 relaxation rate. 
 86      @type R1:                   numpy float array of rank [NE][NS][NM][NO][ND] 
 87      @keyword spin_lock_fields2: The R1rho spin-lock field strengths squared (in rad^2.s^-2). 
 88      @type spin_lock_fields2:    numpy float array of rank [NE][NS][NM][NO][ND] 
 89      @keyword back_calc:         The array for holding the back calculated R1rho values.  Each element corresponds to the combination of theta and spin lock field. 
 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 numer is zero. 
 94      t_numer_zero = False 
 95      t_denom_zero = False 
 96   
 97      # The non-Rex factors. 
 98      sin_theta2 = sin(theta)**2 
 99      R1_R2 = R1 * cos(theta)**2 + r1rho_prime * sin_theta2 
100   
101      # The numerator. 
102      numer = sin_theta2 * phi_ex * kex 
103   
104      # Catch zeros (to avoid pointless mathematical operations). 
105      # This will result in no exchange, returning flat lines. 
106      if min(numer) == 0.0: 
107          t_numer_zero = True 
108          mask_numer_zero = masked_where(numer == 0.0, numer) 
109   
110      # Denominator. 
111      denom = kex**2 + spin_lock_fields2 
112   
113      # Catch math domain error of dividing with 0. 
114      # This is when denom =0. 
115      mask_denom_zero = denom == 0.0 
116      if any(mask_denom_zero): 
117          t_denom_zero = True 
118          denom[mask_denom_zero] = 1.0 
119   
120      # R1rho calculation. 
121      back_calc[:] = R1_R2 + numer / denom 
122   
123      # Replace data in array. 
124      # If numer is zero. 
125      if t_numer_zero: 
126          back_calc[mask_numer_zero.mask] = R1_R2[mask_numer_zero.mask] 
127   
128      # If denom is zero. 
129      if t_denom_zero: 
130          back_calc[mask_denom_zero] = 1e100 
131   
132      # Catch errors, taking a sum over array is the fastest way to check for 
133      # +/- inf (infinity) and nan (not a number). 
134      if not isfinite(sum(back_calc)): 
135          # Replaces nan, inf, etc. with fill value. 
136          fix_invalid(back_calc, copy=False, fill_value=1e100) 
137