Source Code for Module lib.dispersion.tap03

  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 Trott, Abergel and Palmer (2003) off-resonance 2-site exchange R1rho U{TAP03<http://wiki.nmr-relax.com/TAP03>} model. 
 27   
 28  Description 
 29  =========== 
 30   
 31  This module is for the function, gradient and Hessian of the U{TAP03<http://wiki.nmr-relax.com/TAP03>} model. 
 32   
 33   
 34  References 
 35  ========== 
 36   
 37  The model is named after the reference: 
 38   
 39      - Trott, O., Abergel, D., and Palmer, A. (2003).  An average-magnetization analysis of R1rho relaxation outside of the fast exchange.  I{Mol. Phys.}, B{101}, 753-763.  (U{DOI: 10.1080/0026897021000054826<http://dx.doi.org/10.1080/0026897021000054826>}). 
 40   
 41   
 42  Code origin 
 43  =========== 
 44   
 45  The code was copied from the U{TP02<http://wiki.nmr-relax.com/TP02>} model (via the U{MP05<http://wiki.nmr-relax.com/MP05>} 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 TAP03 model can be found in the: 
 57   
 58      - U{relax wiki<http://wiki.nmr-relax.com/TAP03>}, 
 59      - U{relax manual<http://www.nmr-relax.com/manual/The_TAP03_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#TAP03>}. 
 61  """ 
 62   
 63  # Python module imports. 
 64  from numpy import any, arctan2, isfinite, min, sin, sqrt, sum 
 65  from numpy.ma import fix_invalid, masked_where 
 66   
 67   
 68 -def r1rho_TAP03(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 Trott, Abergel and Palmer (2003) equation. 
 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 it is zero. 
 97      t_gamma_neg = False 
 98      t_numer_zero = False 
 99   
100      # Parameter conversions. 
101      pB = 1.0 - pA 
102   
103      # Repetitive calculations (to speed up calculations). 
104      kex2 = kex**2 
105   
106      # Larmor frequency [s^-1]. 
107      Wa = omega 
108      Wb = omega + dw 
109   
110      # Pop-averaged Larmor frequency [s^-1]. 
111      W = pA*Wa + pB*Wb 
112   
113      # The numerator. 
114      phi_ex = pA * pB * dw**2 
115      numer = phi_ex * kex 
116   
117      # The factors. 
118      # Offset of spin-lock from A. 
119      da = Wa - offset 
120   
121      # Offset of spin-lock from B. 
122      db = Wb - offset 
123   
124      # Offset of spin-lock from pop-average. 
125      d = W - offset 
126   
127      # The gamma factor. 
128      sigma = pB*da + pA*db 
129      sigma2 = sigma**2 
130   
131      gamma = 1.0 + phi_ex*(sigma2 - kex2 + spin_lock_fields2) / (sigma2 + kex2 + spin_lock_fields2)**2 
132   
133      # Bad gamma. 
134      mask_gamma_neg = gamma < 0.0 
135      if any(gamma): 
136          t_gamma_neg = True 
137          gamma[mask_gamma_neg] = 0.0 
138   
139      # Special omega values. 
140   
141      # Effective field at A. 
142      waeff2 = gamma*spin_lock_fields2 + da**2 
143   
144      # Effective field at B. 
145      wbeff2 = gamma*spin_lock_fields2 + db**2 
146   
147      # Effective field at pop-average. 
148      weff2 = gamma*spin_lock_fields2 + d**2 
149   
150      # The rotating frame flip angle. 
151      theta = arctan2(spin_lock_fields, d) 
152      hat_theta = arctan2(sqrt(gamma)*spin_lock_fields, d) 
153   
154      # Repetitive calculations (to speed up calculations). 
155      sin_theta2 = sin(theta)**2 
156      hat_sin_theta2 = sin(hat_theta)**2 
157      R1_cos_theta2 = R1 * (1.0 - sin_theta2) 
158      R1rho_prime_sin_theta2 = r1rho_prime * sin_theta2 
159   
160      # Catch zeros (to avoid pointless mathematical operations). 
161      # This will result in no exchange, returning flat lines. 
162      if min(numer) == 0.0: 
163          t_numer_zero = True 
164          mask_numer_zero = masked_where(numer == 0.0, numer) 
165   
166      # Denominator. 
167      denom = waeff2*wbeff2/weff2 + kex2 - 2.0*hat_sin_theta2*phi_ex + (1.0 - gamma)*spin_lock_fields2 
168   
169      # R1rho calculation. 
170      back_calc[:] = R1_cos_theta2 + R1rho_prime_sin_theta2 + hat_sin_theta2 * numer / denom / gamma 
171   
172      # Replace data in array. 
173      # If gamma is negative. 
174      if t_gamma_neg: 
175          back_calc[mask_gamma_neg] = 1e100 
176   
177      # If numer is zero. 
178      if t_numer_zero: 
179          replace = R1_cos_theta2 + R1rho_prime_sin_theta2 
180          back_calc[mask_numer_zero.mask] = replace[mask_numer_zero.mask] 
181   
182      # Catch errors, taking a sum over array is the fastest way to check for 
183      # +/- inf (infinity) and nan (not a number). 
184      if not isfinite(sum(back_calc)): 
185          # Replaces nan, inf, etc. with fill value. 
186          fix_invalid(back_calc, copy=False, fill_value=1e100) 
187