Source Code for Module lib.dispersion.lm63_3site

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2009 Sebastien Morin                                          # 
  4  # Copyright (C) 2013-2014 Edward d'Auvergne                                   # 
  5  # Copyright (C) 2014 Troels E. Linnet                                         # 
  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  r"""The Luz and Meiboom (1963) 3-site fast exchange U{LM63 3-site<http://wiki.nmr-relax.com/LM63_3-site>} model. 
 26   
 27  Description 
 28  =========== 
 29   
 30  This module is for the function, gradient and Hessian of the U{LM63 3-site<http://wiki.nmr-relax.com/LM63_3-site>} model. 
 31   
 32   
 33  References 
 34  ========== 
 35   
 36  The model is named after the reference: 
 37   
 38      - Luz, S. and Meiboom S., (1963)  Nuclear Magnetic Resonance study of protolysis of trimethylammonium ion in aqueous solution - order of reaction with respect to solvent, I{J. Chem. Phys}. B{39}, 366-370 (U{DOI: 10.1063/1.1734254<http://dx.doi.org/10.1063/1.1734254>}). 
 39   
 40   
 41  Equations 
 42  ========= 
 43   
 44  The equation used is:: 
 45   
 46                     _3_ 
 47                     \    phi_ex_i   /     4 * nu_cpmg         /     ki      \ \  
 48      R2eff = R20 +   >   -------- * | 1 - -----------  * tanh | ----------- | | . 
 49                     /__     ki      \         ki              \ 4 * nu_cpmg / / 
 50                     i=2 
 51   
 52  For deconvoluting the parameters, see the relax user manual or the reference: 
 53   
 54      - O'Connell, N. E., Grey, M. J., Tang, Y., Kosuri, P., Miloushev, V. Z., Raleigh, D. P., and Palmer, 3rd, A. G. (2009). Partially folded equilibrium intermediate of the villin headpiece HP67 defined by 13C relaxation dispersion. I{J. Biomol. NMR}, B{45}(1-2), 85-98. (U{DOI: 10.1007/s10858-009-9340-0<http://dx.doi.org/10.1007/s10858-009-9340-0>}). 
 55   
 56   
 57  Links 
 58  ===== 
 59   
 60  More information on the LM63 3-site model can be found in the: 
 61   
 62      - U{relax wiki<http://wiki.nmr-relax.com/LM63_3-site>}, 
 63      - U{relax manual<http://www.nmr-relax.com/manual/The_LM63_3_site_fast_exchange_CPMG_model.html>}, 
 64      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#LM63_3-site>}. 
 65  """ 
 66   
 67  # Python module imports. 
 68  from numpy import fabs, min, tanh, isfinite, sum 
 69  from numpy.ma import fix_invalid, masked_where 
 70   
 71   
 72 -def r2eff_LM63_3site(r20=None, phi_ex_B=None, phi_ex_C=None, kB=None, kC=None, cpmg_frqs=None, back_calc=None): 
 73      """Calculate the R2eff values for the LM63 3-site model. 
 74   
 75      See the module docstring for details. 
 76   
 77   
 78      @keyword r20:           The R20 parameter value (R2 with no exchange). 
 79      @type r20:              numpy float array of rank [NS][NM][NO][ND] 
 80      @keyword phi_ex_B:      The fast exchange factor between sites A and B (ppm^2) 
 81      @type phi_ex_B:         numpy float array of rank [NS][NM][NO][ND] 
 82      @keyword phi_ex_C:      The fast exchange factor between sites A and C (ppm^2) 
 83      @type phi_ex_C:         numpy float array of rank [NE][NS][NM][NO][ND] 
 84      @keyword kB:            Approximate chemical exchange rate constant between sites A and B (the exchange rate in rad/s). 
 85      @type kB:               float 
 86      @keyword kC:            Approximate chemical exchange rate constant between sites A and C (the exchange rate in rad/s). 
 87      @type kC:               float 
 88      @keyword cpmg_frqs:     The CPMG nu1 frequencies. 
 89      @type cpmg_frqs:        numpy float array of rank [NE][NS][NM][NO][ND] 
 90      @keyword back_calc:     The array for holding the back calculated R2eff values.  Each element corresponds to one of the CPMG nu1 frequencies. 
 91      @type back_calc:        numpy float array of rank [NE][NS][NM][NO][ND] 
 92      """ 
 93   
 94      # Once off parameter conversions. 
 95      # The phi_ex_B / kB parameter value. 
 96      rex_B = phi_ex_B / kB 
 97   
 98      # The phi_ex_C / kC parameter value. 
 99      rex_C = phi_ex_C / kC 
100   
101      # Approximate chemical exchange rate constant between sites A and B (the exchange rate in rad/s) divided by 4. 
102      quart_kB = kB / 4.0 
103   
104      # Approximate chemical exchange rate constant between sites A and C (the exchange rate in rad/s) divided by 4. 
105      quart_kC = kC / 4.0 
106   
107      # Flag to tell if values should be replaced. 
108      t_rex_zero = False 
109      t_quart_kB_zero = False 
110      t_quart_kC_zero = False 
111      t_quart_kB_kC_zero = False 
112   
113      # Avoid divide by zero. 
114      if quart_kB == 0.0: 
115          t_quart_kB_zero = True 
116   
117      if quart_kC == 0.0: 
118          t_quart_kC_zero = True 
119   
120      # Test it both is zero. 
121      if t_quart_kB_zero and t_quart_kC_zero: 
122          t_quart_kB_kC_zero = True 
123   
124      # Test if rex 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. 
125      if min(fabs(rex_B)) == 0.0 and min(fabs(rex_C)) == 0.0: 
126          t_rex_zero = True 
127          mask_rex_B_zero = masked_where(rex_B == 0.0, rex_B) 
128          mask_rex_C_zero = masked_where(rex_C == 0.0, rex_C) 
129   
130      # Replace data in array. 
131      # If both quart_kB and quart_kC is zero. 
132      if t_quart_kB_kC_zero: 
133          back_calc[:] = r20 
134      elif t_quart_kB_zero: 
135          back_calc[:] = r20 + rex_C * (1.0 - cpmg_frqs * tanh(quart_kC / cpmg_frqs) / quart_kC) 
136      elif t_quart_kC_zero: 
137          back_calc[:] = r20 + rex_B * (1.0 - cpmg_frqs * tanh(quart_kB / cpmg_frqs) / quart_kB) 
138      else: 
139          # Calc R2eff. 
140          back_calc[:] = r20 + rex_B * (1.0 - cpmg_frqs * tanh(quart_kB / cpmg_frqs) / quart_kB) + rex_C * (1.0 - cpmg_frqs * tanh(quart_kC / cpmg_frqs) / quart_kC) 
141   
142      # If rex is zero. 
143      if t_rex_zero: 
144          back_calc[mask_rex_B_zero.mask] = r20[mask_rex_B_zero.mask] 
145          back_calc[mask_rex_C_zero.mask] = r20[mask_rex_C_zero.mask] 
146   
147      # Catch errors, taking a sum over array is the fastest way to check for 
148      # +/- inf (infinity) and nan (not a number). 
149      if not isfinite(sum(back_calc)): 
150          # Replaces nan, inf, etc. with fill value. 
151          fix_invalid(back_calc, copy=False, fill_value=1e100) 
152