Source Code for Module lib.dispersion.mp05

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2000-2001 Nikolai Skrynnikov                                  # 
  4  # Copyright (C) 2000-2001 Martin Tollinger                                    # 
  5  # Copyright (C) 2013-2014 Edward d'Auvergne                                   # 
  6  # Copyright (C) 2014 Troels E. Linnet                                         # 
  7  #                                                                             # 
  8  # This file is part of the program relax (http://www.nmr-relax.com).          # 
  9  #                                                                             # 
 10  # This program is free software: you can redistribute it and/or modify        # 
 11  # it under the terms of the GNU General Public License as published by        # 
 12  # the Free Software Foundation, either version 3 of the License, or           # 
 13  # (at your option) any later version.                                         # 
 14  #                                                                             # 
 15  # This program is distributed in the hope that it will be useful,             # 
 16  # but WITHOUT ANY WARRANTY without even the implied warranty of               # 
 17  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 18  # GNU General Public License for more details.                                # 
 19  #                                                                             # 
 20  # You should have received a copy of the GNU General Public License           # 
 21  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 22  #                                                                             # 
 23  ############################################################################### 
 24   
 25  # Module docstring. 
 26  """The Miloushev and Palmer (2005) 2-site exchange R1rho U{MP05<http://wiki.nmr-relax.com/MP05>} model. 
 27   
 28  Description 
 29  =========== 
 30   
 31  This module is for the function, gradient and Hessian of the U{MP05<http://wiki.nmr-relax.com/MP05>} model. 
 32   
 33   
 34  References 
 35  ========== 
 36   
 37  The model is named after the reference: 
 38   
 39      - Miloushev V. Z. and Palmer A. (2005).  R1rho relaxation for two-site chemical exchange:  General approximations and some exact solutions.  J. Magn. Reson., 177, 221-227.  (U{DOI: 10.1016/j.jmr.2005.07.023<http://dx.doi.org/10.1016/j.jmr.2005.07.023>}). 
 40   
 41   
 42  Code origin 
 43  =========== 
 44   
 45  The code was copied from the U{TP02<http://wiki.nmr-relax.com/TP02>} model, hence it originates as the funTrottPalmer.m Matlab file from the sim_all.tar file attached to task #7712 (U{https://web.archive.org/web/https://gna.org/task/?7712}).  This is code from Nikolai Skrynnikov and Martin Tollinger. 
 46   
 47  Links to the copyright licensing agreements from all authors are: 
 48   
 49      - Nikolai Skrynnikov, U{http://article.gmane.org/gmane.science.nmr.relax.devel/4279}, 
 50      - Martin Tollinger, U{http://article.gmane.org/gmane.science.nmr.relax.devel/4276}. 
 51   
 52   
 53  Links 
 54  ===== 
 55   
 56  More information on the MP05 model can be found in the: 
 57   
 58      - U{relax wiki<http://wiki.nmr-relax.com/MP05>}, 
 59      - U{relax manual<http://www.nmr-relax.com/manual/The_MP05_2_site_exchange_R1_rho_model.html>}, 
 60      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#MP05>}. 
 61  """ 
 62   
 63  # Python module imports. 
 64  from numpy import arctan2, isfinite, min, sin, sum 
 65  from numpy.ma import fix_invalid, masked_where 
 66   
 67   
 68 -def r1rho_MP05(r1rho_prime=None, omega=None, offset=None, pA=None, dw=None, kex=None, R1=0.0, spin_lock_fields=None, spin_lock_fields2=None, back_calc=None): 
 69      """Calculate the R1rho' values for the TP02 model. 
 70   
 71      See the module docstring for details.  This is the Miloushev and Palmer (2005) equation according to Korzhnev (J. Biomol. NMR (2003), 26, 39-48). 
 72   
 73   
 74      @keyword r1rho_prime:       The R1rho_prime parameter value (R1rho with no exchange). 
 75      @type r1rho_prime:          numpy float array of rank [NE][NS][NM][NO][ND] 
 76      @keyword omega:             The chemical shift for the spin in rad/s. 
 77      @type omega:                numpy float array of rank [NE][NS][NM][NO][ND] 
 78      @keyword offset:            The spin-lock offsets for the data. 
 79      @type offset:               numpy float array of rank [NE][NS][NM][NO][ND] 
 80      @keyword pA:                The population of state A. 
 81      @type pA:                   float 
 82      @keyword dw:                The chemical exchange difference between states A and B in rad/s. 
 83      @type dw:                   numpy float array of rank [NE][NS][NM][NO][ND] 
 84      @keyword kex:               The kex parameter value (the exchange rate in rad/s). 
 85      @type kex:                  float 
 86      @keyword R1:                The R1 relaxation rate. 
 87      @type R1:                   numpy float array of rank [NE][NS][NM][NO][ND] 
 88      @keyword spin_lock_fields:  The R1rho spin-lock field strengths (in rad.s^-1). 
 89      @type spin_lock_fields:     numpy float array of rank [NE][NS][NM][NO][ND] 
 90      @keyword spin_lock_fields2: The R1rho spin-lock field strengths squared (in rad^2.s^-2).  This is for speed. 
 91      @type spin_lock_fields2:    numpy float array of rank [NE][NS][NM][NO][ND] 
 92      @keyword back_calc:         The array for holding the back calculated R1rho values.  Each element corresponds to the combination of offset and spin lock field. 
 93      @type back_calc:            numpy float array of rank [NE][NS][NM][NO][ND] 
 94      """ 
 95   
 96      # Flag to tell if values should be replaced if numer is zero. 
 97      t_numer_zero = False 
 98   
 99      # Parameter conversions. 
100      pB = 1.0 - pA 
101   
102      # Repetitive calculations (to speed up calculations). 
103      kex2 = kex**2 
104   
105      # Larmor frequency [s^-1]. 
106      Wa = omega 
107      Wb = omega + dw 
108   
109      # The numerator. 
110      phi_ex = pA * pB * dw**2 
111      numer = phi_ex * kex 
112   
113      # We assume that A resonates at 0 [s^-1], without loss of generality. 
114      # Pop-averaged Larmor frequency [s^-1]. 
115      W = pA*Wa + pB*Wb 
116   
117      # Offset of spin-lock from A. 
118      da = Wa - offset 
119   
120      # Offset of spin-lock from B. 
121      db = Wb - offset 
122   
123      # Offset of spin-lock from pop-average. 
124      d = W - offset 
125   
126      # Effective field at A. 
127      waeff2 = spin_lock_fields2 + da**2 
128   
129      # Effective field at B. 
130      wbeff2 = spin_lock_fields2 + db**2 
131   
132      # Effective field at pop-average. 
133      weff2 = spin_lock_fields2 + d**2 
134   
135      # The rotating frame flip angle. 
136      theta = arctan2(spin_lock_fields, d) 
137   
138      # Repetitive calculations (to speed up calculations). 
139      sin_theta2 = sin(theta)**2 
140      R1_cos_theta2 = R1 * (1.0 - sin_theta2) 
141      R1rho_prime_sin_theta2 = r1rho_prime * sin_theta2 
142   
143      # Catch zeros (to avoid pointless mathematical operations). 
144      # This will result in no exchange, returning flat lines. 
145      if min(numer) == 0.0: 
146          t_numer_zero = True 
147          mask_numer_zero = masked_where(numer == 0.0, numer) 
148   
149      # Denominator. 
150      waeff2_wbeff2 = waeff2*wbeff2 
151      fact_denom = waeff2_wbeff2 + weff2*kex2 
152   
153      fact = 1.0 + 2.0*kex2*(pA*waeff2 + pB*wbeff2) / fact_denom 
154      denom = waeff2_wbeff2/weff2 + kex2 - sin_theta2*phi_ex*(fact) 
155   
156      # R1rho calculation. 
157      back_calc[:] = R1_cos_theta2 + R1rho_prime_sin_theta2 + sin_theta2 * numer / denom 
158   
159      # Replace data in array. 
160      # If numer is zero. 
161      if t_numer_zero: 
162          replace = R1_cos_theta2 + R1rho_prime_sin_theta2 
163          back_calc[mask_numer_zero.mask] = replace[mask_numer_zero.mask] 
164   
165      # Catch errors, taking a sum over array is the fastest way to check for 
166      # +/- inf (infinity) and nan (not a number). 
167      if not isfinite(sum(back_calc)): 
168          # Replaces nan, inf, etc. with fill value. 
169          fix_invalid(back_calc, copy=False, fill_value=1e100) 
170