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  #                                                                             # 
  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 Trott and Palmer (2002) 2-site exchange R1rho U{TP02<http://wiki.nmr-relax.com/TP02>} model. 
 26   
 27  Description 
 28  =========== 
 29   
 30  This module is for the function, gradient and Hessian of the U{TP02<http://wiki.nmr-relax.com/TP02>} model. 
 31   
 32   
 33  References 
 34  ========== 
 35   
 36  The model is named after the reference: 
 37   
 38      - 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>}). 
 39   
 40   
 41  Code origin 
 42  =========== 
 43   
 44  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. 
 45   
 46  Links to the copyright licensing agreements from all authors are: 
 47   
 48      - Nikolai Skrynnikov, U{http://article.gmane.org/gmane.science.nmr.relax.devel/4279}, 
 49      - Martin Tollinger, U{http://article.gmane.org/gmane.science.nmr.relax.devel/4276}. 
 50   
 51   
 52  Links 
 53  ===== 
 54   
 55  More information on the TP02 model can be found in the: 
 56   
 57      - U{relax wiki<http://wiki.nmr-relax.com/TP02>}, 
 58      - U{relax manual<http://www.nmr-relax.com/manual/TP02_2_site_exchange_R1_model.html>}, 
 59      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#TP02>}. 
 60  """ 
 61   
 62  # Python module imports. 
 63  from math import atan, cos, pi, sin 
 64   
 65   
 66 -def r1rho_TP02(r1rho_prime=None, omega=None, offset=None, pA=None, pB=None, dw=None, kex=None, R1=0.0, spin_lock_fields=None, spin_lock_fields2=None, back_calc=None, num_points=None): 
 67      """Calculate the R1rho' values for the TP02 model. 
 68   
 69      See the module docstring for details.  This is the Trott and Palmer (2002) equation according to Korzhnev (J. Biomol. NMR (2003), 26, 39-48). 
 70   
 71   
 72      @keyword r1rho_prime:       The R1rho_prime parameter value (R1rho with no exchange). 
 73      @type r1rho_prime:          float 
 74      @keyword omega:             The chemical shift for the spin in rad/s. 
 75      @type omega:                float 
 76      @keyword offset:            The spin-lock offsets for the data. 
 77      @type offset:               numpy rank-1 float array 
 78      @keyword pA:                The population of state A. 
 79      @type pA:                   float 
 80      @keyword pB:                The population of state B. 
 81      @type pB:                   float 
 82      @keyword dw:                The chemical exchange difference between states A and B in rad/s. 
 83      @type dw:                   float 
 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:                   float 
 88      @keyword spin_lock_fields:  The R1rho spin-lock field strengths (in rad.s^-1). 
 89      @type spin_lock_fields:     numpy rank-1 float array 
 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 rank-1 float array 
 92      @keyword back_calc:         The array for holding the back calculated R1rho values.  Each element corresponds to one of the spin-lock fields. 
 93      @type back_calc:            numpy rank-1 float array 
 94      @keyword num_points:        The number of points on the dispersion curve, equal to the length of the spin_lock_fields and back_calc arguments. 
 95      @type num_points:           int 
 96      """ 
 97   
 98      # Repetitive calculations (to speed up calculations). 
 99      Wa = omega                  # Larmor frequency [s^-1]. 
100      Wb = omega + dw             # Larmor frequency [s^-1]. 
101      kex2 = kex**2 
102      W = pA*Wa + pB*Wb           # Pop-averaged Larmor frequency [s^-1]. 
103      da = Wa - offset            # Offset of spin-lock from A. 
104      db = Wb - offset            # Offset of spin-lock from B. 
105      d = W - offset              # Offset of spin-lock from pop-average. 
106      da2 = da**2 
107      db2 = db**2 
108      d2 = d**2 
109   
110      # The numerator. 
111      numer = pA * pB * dw**2 * kex 
112   
113      # Loop over the dispersion points, back calculating the R1rho values. 
114      for i in range(num_points): 
115          # We assume that A resonates at 0 [s^-1], without loss of generality. 
116          waeff2 = spin_lock_fields2[i] + da2       # Effective field at A. 
117          wbeff2 = spin_lock_fields2[i] + db2       # Effective field at B. 
118          weff2 = spin_lock_fields2[i] + d2         # Effective field at pop-average. 
119   
120          # The rotating frame flip angle. 
121          theta = atan(spin_lock_fields[i] / d) 
122   
123          # Repetitive calculations (to speed up calculations). 
124          sin_theta2 = sin(theta)**2 
125          R1_cos_theta2 = R1 * (1.0 - sin_theta2) 
126          R1rho_prime_sin_theta2 = r1rho_prime * sin_theta2 
127   
128          # Catch zeros (to avoid pointless mathematical operations). 
129          if numer == 0.0: 
130              back_calc[i] = R1_cos_theta2 + R1rho_prime_sin_theta2 
131              continue 
132   
133          # Denominator. 
134          denom = waeff2 * wbeff2 / weff2 + kex2 
135          #denom_extended = waeff2*wbeff2/weff2+kex2-2*sin_theta2*pA*pB*dw**2 
136    
137          # Avoid divide by zero. 
138          if denom == 0.0: 
139              back_calc[i] = 1e100 
140              continue 
141   
142          # R1rho calculation. 
143          back_calc[i] = R1_cos_theta2 + R1rho_prime_sin_theta2 + sin_theta2 * numer / denom 
144