Source Code for Module lib.dispersion.m61b

  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 Meiboom (1961) 2-site on-resonance skewed population R1rho U{M61 skew<http://wiki.nmr-relax.com/M61_skew>} model. 
 25   
 26  Description 
 27  =========== 
 28   
 29  This module is for the function, gradient and Hessian of the U{M61 skew<http://wiki.nmr-relax.com/M61_skew>} model. 
 30   
 31   
 32  References 
 33  ========== 
 34   
 35  The model is named after the reference: 
 36   
 37      - Meiboom S. (1961).  Nuclear magnetic resonance study of the proton transfer in water.  I{J. Chem. Phys.}, B{34}, 375-388.  (U{DOI: 10.1063/1.1700960<http://dx.doi.org/10.1063/1.1700960>}). 
 38   
 39   
 40  Equations 
 41  ========= 
 42   
 43  The equation used is:: 
 44   
 45                             pA^2.pB.delta_omega^2.kex 
 46      R1rho = R1rho' + -------------------------------------- , 
 47                       kex^2 + pA^2.delta_omega^2 + omega_1^2 
 48   
 49  where R1rho' is the R1rho value in the absence of exchange, 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_1 = omega_e is the effective field in the rotating frame. 
 50   
 51   
 52  Links 
 53  ===== 
 54   
 55  More information on the M61 skew model can be found in the: 
 56   
 57      - U{relax wiki<http://wiki.nmr-relax.com/M61_skew>}, 
 58      - U{relax manual<http://www.nmr-relax.com/manual/M61_skew_2_site_fast_exchange_R1_model.html>}, 
 59      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#M61_skew>}. 
 60  """ 
 61   
 62  # Python module imports. 
 63  from numpy import abs, array, isfinite, min, sum 
 64   
 65 -def r1rho_M61b(r1rho_prime=None, pA=None, dw=None, kex=None, spin_lock_fields2=None, back_calc=None, num_points=None): 
 66      """Calculate the R1rho values for the M61 skew model. 
 67   
 68      See the module docstring for details. 
 69   
 70   
 71      @keyword r1rho_prime:       The R1rho_prime parameter value (R1rho with no exchange). 
 72      @type r1rho_prime:          float 
 73      @keyword pA:                The population of state A. 
 74      @type pA:                   float 
 75      @keyword dw:                The chemical exchange difference between states A and B in rad/s. 
 76      @type dw:                   float 
 77      @keyword kex:               The kex parameter value (the exchange rate in rad/s). 
 78      @type kex:                  float 
 79      @keyword spin_lock_fields2: The R1rho spin-lock field strengths squared (in rad^2.s^-2). 
 80      @type spin_lock_fields2:    numpy rank-1 float array 
 81      @keyword back_calc:         The array for holding the back calculated R1rho values.  Each element corresponds to the combination of spin lock field. 
 82      @type back_calc:            numpy rank-1 float array 
 83      @keyword num_points:        The number of points on the dispersion curve, equal to the length of the spin_lock_fields and back_calc arguments. 
 84      @type num_points:           int 
 85      """ 
 86   
 87      # The B population. 
 88      pB = 1.0 - pA 
 89   
 90      # Repetitive calculations (to speed up calculations). 
 91      pA2dw2 = pA**2 * dw**2 
 92      kex2_pA2dw2 = kex**2 + pA2dw2 
 93   
 94      # The numerator. 
 95      numer = pA2dw2 * pB * kex 
 96   
 97      # Catch zeros (to avoid pointless mathematical operations). 
 98      # This will result in no exchange, returning flat lines. 
 99      if numer == 0.0: 
100          back_calc[:] = array([r1rho_prime]*num_points) 
101          return 
102   
103      # Denominator. 
104      denom = kex2_pA2dw2 + spin_lock_fields2 
105   
106      # Catch math domain error of dividing with 0. 
107      # This is when denom=0. 
108      if min(abs(denom)) == 0: 
109          back_calc[:] = array([1e100]*num_points) 
110          return 
111   
112   
113      # R1rho calculation. 
114      R1rho = r1rho_prime + numer / denom 
115   
116      # Catch errors, taking a sum over array is the fastest way to check for 
117      # +/- inf (infinity) and nan (not a number). 
118      if not isfinite(sum(R1rho)): 
119          R1rho = array([1e100]*num_points) 
120   
121      back_calc[:] = R1rho 
122