Source Code for Module lib.dispersion.tp02

  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 and Palmer (2002) 2-site exchange R1rho U{TP02<http://wiki.nmr-relax.com/TP02>} model. 
 27   
 28  Description 
 29  =========== 
 30   
 31  This module is for the function, gradient and Hessian of the U{TP02<http://wiki.nmr-relax.com/TP02>} model. 
 32   
 33   
 34  References 
 35  ========== 
 36   
 37  The model is named after the reference: 
 38   
 39      - Trott, O. and Palmer, 3rd, A. G. (2002). R1rho relaxation outside of the fast-exchange limit. I{J. Magn. Reson.}, B{154}(1), 157-160. (U{DOI: 10.1006/jmre.2001.2466<http://dx.doi.org/10.1006/jmre.2001.2466>}). 
 40   
 41   
 42  Code origin 
 43  =========== 
 44   
 45  The code 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 TP02 model can be found in the: 
 57   
 58      - U{relax wiki<http://wiki.nmr-relax.com/TP02>}, 
 59      - U{relax manual<http://www.nmr-relax.com/manual/The_TP02_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#TP02>}. 
 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_TP02(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 and Palmer (2002) 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 it 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      # Pop-averaged Larmor frequency [s^-1]. 
110      W = pA*Wa + pB*Wb 
111   
112      # Offset of spin-lock from A. 
113      da = Wa - offset 
114   
115      # Offset of spin-lock from B. 
116      db = Wb - offset 
117   
118      # Offset of spin-lock from pop-average. 
119      d = W - offset 
120      da2 = da**2 
121      db2 = db**2 
122      d2 = d**2 
123   
124      # The numerator. 
125      numer = pA * pB * dw**2 * kex 
126   
127      # We assume that A resonates at 0 [s^-1], without loss of generality. 
128      # Effective field at A. 
129      waeff2 = spin_lock_fields2 + da2 
130   
131      # Effective field at B. 
132      wbeff2 = spin_lock_fields2 + db2 
133   
134      # Effective field at pop-average. 
135      weff2 = spin_lock_fields2 + d2 
136   
137      # The rotating frame flip angle. 
138      theta = arctan2(spin_lock_fields, d) 
139   
140      # Repetitive calculations (to speed up calculations). 
141      sin_theta2 = sin(theta)**2 
142      R1_cos_theta2 = R1 * (1.0 - sin_theta2) 
143      R1rho_prime_sin_theta2 = r1rho_prime * sin_theta2 
144   
145      # Catch zeros (to avoid pointless mathematical operations). 
146      # This will result in no exchange, returning flat lines. 
147      if min(numer) == 0.0: 
148          t_numer_zero = True 
149          mask_numer_zero = masked_where(numer == 0.0, numer) 
150   
151      # Denominator. 
152      denom = waeff2 * wbeff2 / weff2 + kex2 
153      #denom_extended = waeff2*wbeff2/weff2+kex2-2*sin_theta2*pA*pB*dw**2 
154   
155      # R1rho calculation. 
156      back_calc[:] = R1_cos_theta2 + R1rho_prime_sin_theta2 + sin_theta2 * numer / denom 
157   
158      # If numer is zero. 
159      if t_numer_zero: 
160          replace = R1_cos_theta2 + R1rho_prime_sin_theta2 
161          back_calc[mask_numer_zero.mask] = replace[mask_numer_zero.mask] 
162   
163      # Catch errors, taking a sum over array is the fastest way to check for 
164      # +/- inf (infinity) and nan (not a number). 
165      if not isfinite(sum(back_calc)): 
166          # Replaces nan, inf, etc. with fill value. 
167          fix_invalid(back_calc, copy=False, fill_value=1e100) 
168