Package lib :: Package dispersion :: Module m61
[hide private]
[frames] | no frames]

Source Code for Module lib.dispersion.m61

  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  """The Meiboom (1961) 2-site fast exchange R1rho U{M61<http://wiki.nmr-relax.com/M61>} model. 
 26   
 27  Description 
 28  =========== 
 29   
 30  This module is for the function, gradient and Hessian of the U{M61<http://wiki.nmr-relax.com/M61>} model. 
 31   
 32   
 33  References 
 34  ========== 
 35   
 36  The model is named after the reference: 
 37   
 38      - Meiboom S. (1961).  Nuclear magnetic resonance study of the proton transfer in water.  I{J. Chem. Phys.}, B{34}, 375-388.  (U{DOI: 10.1063/1.1700960<http://dx.doi.org/10.1063/1.1700960>}). 
 39   
 40   
 41  Equations 
 42  ========= 
 43   
 44  The equation used is:: 
 45   
 46                         phi_ex * kex 
 47      R1rho = R1rho' + ----------------- , 
 48                       kex^2 + omega_1^2 
 49   
 50  where:: 
 51   
 52      phi_ex = pA * pB * delta_omega^2 , 
 53   
 54  R1rho' is the R1rho value in the absence of exchange, kex is the chemical exchange rate constant, pA and pB are the populations of states A and B, delta_omega is the chemical shift difference between the two states, and omega_1 is the spin-lock field strength. 
 55   
 56   
 57  Links 
 58  ===== 
 59   
 60  More information on the M61 model can be found in the: 
 61   
 62      - U{relax wiki<http://wiki.nmr-relax.com/M61>}, 
 63      - U{relax manual<http://www.nmr-relax.com/manual/The_M61_2_site_fast_exchange_R1_rho_model.html>}, 
 64      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#M61>}. 
 65  """ 
 66   
 67  # Python module imports. 
 68  from numpy import any, isfinite, min, sum 
 69  from numpy.ma import fix_invalid, masked_where 
 70   
 71   
72 -def r1rho_M61(r1rho_prime=None, phi_ex=None, kex=None, spin_lock_fields2=None, back_calc=None):
73 """Calculate the R2eff values for the M61 model. 74 75 See the module docstring for details. 76 77 78 @keyword r1rho_prime: The R1rho_prime parameter value (R1rho with no exchange). 79 @type r1rho_prime: numpy float array of rank [NE][NS][NM][NO][ND] 80 @keyword phi_ex: The phi_ex parameter value (pA * pB * delta_omega^2). 81 @type phi_ex: numpy float array of rank [NE][NS][NM][NO][ND] 82 @keyword kex: The kex parameter value (the exchange rate in rad/s). 83 @type kex: float 84 @keyword spin_lock_fields2: The R1rho spin-lock field strengths squared (in rad^2.s^-2). 85 @type spin_lock_fields2: numpy float array of rank [NE][NS][NM][NO][ND] 86 @keyword back_calc: The array for holding the back calculated R1rho values. Each element corresponds to the combination of spin lock field. 87 @type back_calc: numpy float array of rank [NE][NS][NM][NO][ND] 88 """ 89 90 # Flag to tell if values should be replaced if numer is zero. 91 t_numer_zero = False 92 t_denom_zero = False 93 94 # Repetitive calculations (to speed up calculations). 95 kex2 = kex**2 96 97 # The numerator. 98 numer = phi_ex * kex 99 100 # Catch zeros (to avoid pointless mathematical operations). 101 # This will result in no exchange, returning flat lines. 102 if min(numer) == 0.0: 103 t_numer_zero = True 104 mask_numer_zero = masked_where(numer == 0.0, numer) 105 106 # Denominator. 107 denom = kex2 + spin_lock_fields2 108 109 # Catch math domain error of dividing with 0. 110 # This is when denom = 0. 111 mask_denom_zero = denom == 0.0 112 if any(mask_denom_zero): 113 t_denom_zero = True 114 denom[mask_denom_zero] = 1.0 115 116 # R1rho calculation. 117 back_calc[:] = r1rho_prime + numer / denom 118 119 # Replace data in array. 120 # If numer is zero. 121 if t_numer_zero: 122 back_calc[mask_numer_zero.mask] = r1rho_prime[mask_numer_zero.mask] 123 124 # If denom is zero. 125 if t_denom_zero: 126 back_calc[mask_denom_zero] = 1e100 127 128 # Catch errors, taking a sum over array is the fastest way to check for 129 # +/- inf (infinity) and nan (not a number). 130 if not isfinite(sum(back_calc)): 131 # Replaces nan, inf, etc. with fill value. 132 fix_invalid(back_calc, copy=False, fill_value=1e100)
133