Source Code for Module lib.dispersion.mmq_cr72

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2013-2014 Edward d'Auvergne                                   # 
  4  #                                                                             # 
  5  # This file is part of the program relax (http://www.nmr-relax.com).          # 
  6  #                                                                             # 
  7  # This program is free software: you can redistribute it and/or modify        # 
  8  # it under the terms of the GNU General Public License as published by        # 
  9  # the Free Software Foundation, either version 3 of the License, or           # 
 10  # (at your option) any later version.                                         # 
 11  #                                                                             # 
 12  # This program is distributed in the hope that it will be useful,             # 
 13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 15  # GNU General Public License for more details.                                # 
 16  #                                                                             # 
 17  # You should have received a copy of the GNU General Public License           # 
 18  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 19  #                                                                             # 
 20  ############################################################################### 
 21   
 22  # Module docstring. 
 23  """The CR72 model extended for MMQ CPMG data, called the U{MMQ CR72<http://wiki.nmr-relax.com/MMQ_CR72>} model. 
 24   
 25  Description 
 26  =========== 
 27   
 28  This module is for the function, gradient and Hessian of the U{MMQ CR72<http://wiki.nmr-relax.com/MMQ_CR72>} model. 
 29   
 30   
 31  References 
 32  ========== 
 33   
 34  The Carver and Richards (1972) 2-site model for all times scales was extended for multiple-MQ (MMQ) CPMG data by: 
 35   
 36      - Korzhnev, D. M., Kloiber, K., Kanelis, V., Tugarinov, V., and Kay, L. E. (2004).  Probing slow dynamics in high molecular weight proteins by methyl-TROSY NMR spectroscopy: Application to a 723-residue enzyme.  I{J. Am. Chem. Soc.}, B{126}, 3964-3973.  (U{DOI: 10.1021/ja039587i<http://dx.doi.org/10.1021/ja039587i>}). 
 37   
 38   
 39  Links 
 40  ===== 
 41   
 42  More information on the MMQ CR72 model can be found in the: 
 43   
 44      - U{relax wiki<http://wiki.nmr-relax.com/MMQ_CR72>}, 
 45      - U{relax manual<http://www.nmr-relax.com/manual/MMQ_CR72_model.html>}, 
 46      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#MMQ_CR72>}. 
 47  """ 
 48   
 49  # Python module imports. 
 50  from numpy import arccosh, array, cos, cosh, isfinite, log, max, sin, sqrt, sum 
 51   
 52   
 53 -def r2eff_mmq_cr72(r20=None, pA=None, pB=None, dw=None, dwH=None, kex=None, k_AB=None, k_BA=None, cpmg_frqs=None, inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None): 
 54      """The CR72 model extended to MMQ CPMG data. 
 55   
 56      This function calculates and stores the R2eff values. 
 57   
 58   
 59      @keyword r20:           The R2 value in the absence of exchange. 
 60      @type r20:              float 
 61      @keyword pA:            The population of state A. 
 62      @type pA:               float 
 63      @keyword pB:            The population of state B. 
 64      @type pB:               float 
 65      @keyword dw:            The chemical exchange difference between states A and B in rad/s. 
 66      @type dw:               float 
 67      @keyword dwH:           The proton chemical exchange difference between states A and B in rad/s. 
 68      @type dwH:              float 
 69      @keyword kex:           The kex parameter value (the exchange rate in rad/s). 
 70      @type kex:              float 
 71      @keyword k_AB:          The rate of exchange from site A to B (rad/s). 
 72      @type k_AB:             float 
 73      @keyword k_BA:          The rate of exchange from site B to A (rad/s). 
 74      @type k_BA:             float 
 75      @keyword cpmg_frqs:     The CPMG nu1 frequencies. 
 76      @type cpmg_frqs:        numpy rank-1 float array 
 77      @keyword inv_tcpmg:     The inverse of the total duration of the CPMG element (in inverse seconds). 
 78      @type inv_tcpmg:        float 
 79      @keyword tcp:           The tau_CPMG times (1 / 4.nu1). 
 80      @type tcp:              numpy rank-1 float array 
 81      @keyword back_calc:     The array for holding the back calculated R2eff values.  Each element corresponds to one of the CPMG nu1 frequencies. 
 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 tcp and back_calc arguments. 
 84      @type num_points:       int 
 85      @keyword power:         The matrix exponential power array. 
 86      @type power:            numpy int16, rank-1 array 
 87      """ 
 88   
 89      # Catch parameter values that will result in no exchange, returning flat R2eff = R20 lines (when kex = 0.0, k_AB = 0.0). 
 90      if (dw == 0.0 and dwH == 0.0) or pA == 1.0 or k_AB == 0.0: 
 91          back_calc[:] = array([r20]*num_points) 
 92          return 
 93   
 94      # Repetitive calculations (to speed up calculations). 
 95      dw2 = dw**2 
 96      r20_kex = r20 + kex/2.0 
 97      pApBkex2 = k_AB * k_BA 
 98      isqrt_pApBkex2 = 1.j*sqrt(pApBkex2) 
 99      sqrt_pBpA = sqrt(pB/pA) 
100      ikex = 1.j*kex 
101   
102      # The d+/- values. 
103      d = dwH + dw 
104      dpos = d + ikex 
105      dneg = d - ikex 
106   
107      # The z+/- values. 
108      z = dwH - dw 
109      zpos = z + ikex 
110      zneg = z - ikex 
111   
112      # The Psi and zeta values. 
113      fact = 1.j*dwH + k_BA - k_AB 
114      Psi = fact**2 - dw2 + 4.0*pApBkex2 
115      zeta = -2.0*dw * fact 
116   
117      # More repetitive calculations. 
118      sqrt_psi2_zeta2 = sqrt(Psi**2 + zeta**2) 
119   
120      # The D+/- values. 
121      D_part = (Psi + 2.0*dw2) / sqrt_psi2_zeta2 
122      Dpos = 0.5 * (1.0 + D_part) 
123      Dneg = 0.5 * (-1.0 + D_part) 
124   
125      # Partial eta+/- values. 
126      eta_scale = 2.0**(-3.0/2.0) 
127      etapos_part = eta_scale * sqrt(Psi + sqrt_psi2_zeta2) 
128      etaneg_part = eta_scale * sqrt(-Psi + sqrt_psi2_zeta2) 
129   
130      # The full eta+ values. 
131      etapos = etapos_part / cpmg_frqs 
132   
133      # Catch math domain error of cosh(val > 710). 
134      # This is when etapos > 710. 
135      if max(etapos) > 700: 
136          back_calc[:] = array([r20]*num_points) 
137          return 
138   
139      # The full eta - values. 
140      etaneg = etaneg_part / cpmg_frqs 
141   
142      # The mD value. 
143      mD = isqrt_pApBkex2 / (dpos * zpos) * (zpos + 2.0*dw*sin(zpos*tcp)/sin((dpos + zpos)*tcp)) 
144   
145      # The mZ value. 
146      mZ = -isqrt_pApBkex2 / (dneg * zneg) * (dneg - 2.0*dw*sin(dneg*tcp)/sin((dneg + zneg)*tcp)) 
147   
148      # The Q value. 
149      Q = 1 - mD**2 + mD*mZ - mZ**2 + 0.5*(mD + mZ)*sqrt_pBpA 
150      Q = Q.real 
151   
152      # The first eigenvalue. 
153      lambda1 = r20_kex - cpmg_frqs * arccosh(Dpos * cosh(etapos) - Dneg * cos(etaneg)) 
154   
155      # The full formula. 
156      R2eff = lambda1.real - inv_tcpmg * log(Q) 
157   
158      # Catch errors, taking a sum over array is the fastest way to check for 
159      # +/- inf (infinity) and nan (not a number). 
160      if not isfinite(sum(R2eff)): 
161          R2eff = array([1e100]*num_points) 
162   
163      back_calc[:] = R2eff 
164