Source Code for Module lib.dispersion.mmq_cr72

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2013-2014 Edward d'Auvergne                                   # 
  4  # Copyright (C) 2014 Troels E. Linnet                                         # 
  5  #                                                                             # 
  6  # This file is part of the program relax (http://www.nmr-relax.com).          # 
  7  #                                                                             # 
  8  # This program is free software: you can redistribute it and/or modify        # 
  9  # it under the terms of the GNU General Public License as published by        # 
 10  # the Free Software Foundation, either version 3 of the License, or           # 
 11  # (at your option) any later version.                                         # 
 12  #                                                                             # 
 13  # This program is distributed in the hope that it will be useful,             # 
 14  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 15  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 16  # GNU General Public License for more details.                                # 
 17  #                                                                             # 
 18  # You should have received a copy of the GNU General Public License           # 
 19  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 20  #                                                                             # 
 21  ############################################################################### 
 22   
 23  # Module docstring. 
 24  """The CR72 model extended for MMQ CPMG data, called the U{MMQ CR72<http://wiki.nmr-relax.com/MMQ_CR72>} model. 
 25   
 26  Description 
 27  =========== 
 28   
 29  This module is for the function, gradient and Hessian of the U{MMQ CR72<http://wiki.nmr-relax.com/MMQ_CR72>} model. 
 30   
 31   
 32  References 
 33  ========== 
 34   
 35  The Carver and Richards (1972) 2-site model for all times scales was extended for multiple-MQ (MMQ) CPMG data by: 
 36   
 37      - 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>}). 
 38   
 39   
 40  Links 
 41  ===== 
 42   
 43  More information on the MMQ CR72 model can be found in the: 
 44   
 45      - U{relax wiki<http://wiki.nmr-relax.com/MMQ_CR72>}, 
 46      - U{relax manual<http://www.nmr-relax.com/manual/The_MMQ_CR72_model.html>}, 
 47      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#MMQ_CR72>}. 
 48  """ 
 49   
 50  # Python module imports. 
 51  from numpy import arccosh, cos, cosh, isfinite, fabs, log, min, max, sin, sqrt, sum 
 52  from numpy.ma import fix_invalid, masked_greater_equal, masked_where 
 53   
 54  # Repetitive calculations (to speed up calculations). 
 55  eta_scale = 2.0**(-3.0/2.0) 
 56   
 57   
 58 -def r2eff_mmq_cr72(r20=None, pA=None, dw=None, dwH=None, kex=None, cpmg_frqs=None, inv_tcpmg=None, tcp=None, back_calc=None): 
 59      """The CR72 model extended to MMQ CPMG data. 
 60   
 61      This function calculates and stores the R2eff values. 
 62   
 63   
 64      @keyword r20:           The R2 value in the absence of exchange. 
 65      @type r20:              numpy float array of rank [NS][NM][NO][ND] 
 66      @keyword pA:            The population of state A. 
 67      @type pA:               float 
 68      @keyword dw:            The chemical exchange difference between states A and B in rad/s. 
 69      @type dw:               numpy float array of rank [NS][NM][NO][ND] 
 70      @keyword dwH:           The proton chemical exchange difference between states A and B in rad/s. 
 71      @type dwH:              numpy float array of rank [NS][NM][NO][ND] 
 72      @keyword kex:           The kex parameter value (the exchange rate in rad/s). 
 73      @type kex:              float 
 74      @keyword cpmg_frqs:     The CPMG nu1 frequencies. 
 75      @type cpmg_frqs:        numpy float array of rank [NS][NM][NO][ND] 
 76      @keyword inv_tcpmg:     The inverse of the total duration of the CPMG element (in inverse seconds). 
 77      @type inv_tcpmg:        numpy float array of rank [NS][NM][NO][ND] 
 78      @keyword tcp:           The tau_CPMG times (1 / 4.nu1). 
 79      @type tcp:              numpy float array of rank [NS][NM][NO][ND] 
 80      @keyword back_calc:     The array for holding the back calculated R2eff values.  Each element corresponds to one of the CPMG nu1 frequencies. 
 81      @type back_calc:        numpy float array of rank [NS][NM][NO][ND] 
 82      """ 
 83   
 84      # Once off parameter conversions. 
 85      pB = 1.0 - pA 
 86      k_BA = pA * kex 
 87      k_AB = pB * kex 
 88   
 89      # Flag to tell if values should be replaced if max_etapos in cosh function is violated. 
 90      t_dw_dw_H_zero = False 
 91      t_max_etapos = False 
 92   
 93      # Test if pA or kex is zero. 
 94      if kex == 0.0 or pA == 1.0: 
 95          back_calc[:] = r20 
 96          return 
 97   
 98      # Test if dw and dwH is zero. Create a mask for the affected spins to replace these with R20 at the end of the calculationWait for replacement, since this is spin specific. 
 99      if min(fabs(dw)) == 0.0 and min(fabs(dwH)) == 0.0: 
100          t_dw_dw_H_zero = True 
101          mask_dw_zero = masked_where(dw == 0.0, dw) 
102          mask_dw_H_zero = masked_where(dwH == 0.0, dwH) 
103   
104      # Repetitive calculations (to speed up calculations). 
105      dw2 = dw**2 
106      r20_kex = r20 + kex/2.0 
107      pApBkex2 = k_AB * k_BA 
108      isqrt_pApBkex2 = 1.j*sqrt(pApBkex2) 
109      sqrt_pBpA = sqrt(pB/pA) 
110      ikex = 1.j*kex 
111   
112      # The d+/- values. 
113      d = dwH + dw 
114      dpos = d + ikex 
115      dneg = d - ikex 
116   
117      # The z+/- values. 
118      z = dwH - dw 
119      zpos = z + ikex 
120      zneg = z - ikex 
121   
122      # The Psi and zeta values. 
123      fact = 1.j*dwH + k_BA - k_AB 
124      Psi = fact**2 - dw2 + 4.0*pApBkex2 
125      zeta = -2.0*dw * fact 
126   
127      # More repetitive calculations. 
128      sqrt_psi2_zeta2 = sqrt(Psi**2 + zeta**2) 
129   
130      # The D+/- values. 
131      D_part = (0.5*Psi + dw2) / sqrt_psi2_zeta2 
132      Dpos = 0.5 + D_part 
133      Dneg = -0.5 + D_part 
134   
135      # The eta+/- values. 
136      eta_fact = eta_scale / cpmg_frqs 
137      etapos = eta_fact * sqrt(Psi + sqrt_psi2_zeta2) 
138      etaneg = eta_fact * sqrt(-Psi + sqrt_psi2_zeta2) 
139   
140      # Catch math domain error of cosh(val > 710). 
141      # This is when etapos > 710. 
142      if max(etapos) > 700: 
143          t_max_etapos = True 
144          mask_max_etapos = masked_greater_equal(etapos, 700.0) 
145          # To prevent math errors, set etapos to 1. 
146          etapos[mask_max_etapos.mask] = 1.0 
147   
148      # The mD value. 
149      mD = isqrt_pApBkex2 / (dpos * zpos) * (zpos + 2.0*dw*sin(zpos*tcp)/sin((dpos + zpos)*tcp)) 
150   
151      # The mZ value. 
152      mZ = -isqrt_pApBkex2 / (dneg * zneg) * (dneg - 2.0*dw*sin(dneg*tcp)/sin((dneg + zneg)*tcp)) 
153   
154      # The Q value. 
155      Q = 1 - mD**2 + mD*mZ - mZ**2 + 0.5*(mD + mZ)*sqrt_pBpA 
156      Q = Q.real 
157   
158      # The first eigenvalue. 
159      lambda1 = r20_kex - cpmg_frqs * arccosh(Dpos * cosh(etapos) - Dneg * cos(etaneg)) 
160   
161      # The full formula. 
162      back_calc[:] = lambda1.real - inv_tcpmg * log(Q) 
163   
164      # Replace data in array. 
165      # If eta_pos above 700. 
166      if t_max_etapos: 
167          back_calc[mask_max_etapos.mask] = r20[mask_max_etapos.mask] 
168   
169      # Replace data in array. 
170      # If dw and dwH is zero. 
171      if t_dw_dw_H_zero: 
172          back_calc[mask_dw_zero.mask] = r20[mask_dw_zero.mask] 
173          back_calc[mask_dw_H_zero.mask] = r20[mask_dw_H_zero.mask] 
174   
175      # Catch errors, taking a sum over array is the fastest way to check for 
176      # +/- inf (infinity) and nan (not a number). 
177      if not isfinite(sum(back_calc)): 
178          # Replaces nan, inf, etc. with fill value. 
179          fix_invalid(back_calc, copy=False, fill_value=1e100) 
180