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

Source Code for Module lib.dispersion.b14

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2014 Troels E. Linnet                                         # 
  4  # Copyright (C) 2014 Andrew Baldwin                                           # 
  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 Baldwin (2014) 2-site exact solution model for all time scales U{B14<http://wiki.nmr-relax.com/B14>}. 
 25   
 26  Description 
 27  =========== 
 28   
 29  This module is for the function, gradient and Hessian of the U{B14<http://wiki.nmr-relax.com/B14>} model. 
 30   
 31   
 32  References 
 33  ========== 
 34   
 35  The model is named after the reference: 
 36   
 37      - Andrew J. Baldwin (2014).  An exact solution for R2,eff in CPMG experiments in the case of two site chemical exchange.  I{J. Magn. Reson.}.  (U{DOI: 10.1016/j.jmr.2014.02.023 <http://dx.doi.org/10.1016/j.jmr.2014.02.023>}). 
 38   
 39   
 40  Equations 
 41  ========= 
 42   
 43  The equation used is:: 
 44   
 45              R2A0 + R2B0 + kex      Ncyc                      1      ( 1+y            1-y                          ) 
 46      R2eff = ------------------ -  ------ * cosh^-1 * v1c - ------ ln( --- + ------------------ * (v2 + 2*kAB*pD ) ) , 
 47                    2                Trel                     Trel    (  2    2*sqrt(v1c^2 -1 )                     ) 
 48   
 49                              1      ( 1+y            1-y                          ) 
 50            = R2eff(CR72) - ------ ln( --- + ------------------ * (v2 + 2*kAB*pD ) ) , 
 51                             Trel    (  2    2*sqrt(v1c^2 -1 )                     ) 
 52   
 53  Which have these following definitions:: 
 54   
 55      v1c = F0 * cosh(tauCP * E0)- F2 * cosh(tauCP * E2) , 
 56      v1s = F0 * sinh(tauCP * E0)- F2 * sinh(tauCP * E2) , 
 57      v2*N = v1s * (OB-OA) + 4OB * F1^a * sinh(tauCP * E1) , 
 58      pD N = v1s + (F1^a + F1^b) * sinh(tauCP * E1) , 
 59      v3 = ( v2^2 + 4 * kBA * kAB * pD^2 )^1/2 , 
 60      y = ( (v1c-v3)/(v1c+v3) )^NCYC , 
 61   
 62  Note, E2 is complex. If |x| denotes the complex modulus:: 
 63   
 64      cosh(tauCP * E2) = cos(tauCP * |E2|) , 
 65      sinh(tauCP * E2) = i sin(tauCP * |E2|) , 
 66   
 67  The term pD is based on product of the off diagonal elements in the CPMG propagator (Supplementary Section 3). 
 68   
 69  It is interesting to consider the region of validity of the Carver Richards result.  The two results are equal when the correction is zero, which is true when:: 
 70   
 71      sqrt(v1c^2-1) ~ v2 + 2*kAB * pD , 
 72   
 73  This occurs when 2*kAB * pD tends to zero, and so v2=v3.  Setting "kAB * pD" to zero, amounts to neglecting magnetisation that starts on the ground state ensemble and end on the excited state ensemble and vice versa.  This will be a good approximation when pA >> p_B. 
 74   
 75  In practise, significant deviations from the Carver Richards equation can be incurred if pB > 1 %.  Incorporation of the correction term into equation (50), results in an improved description of the CPMG experiment over the Carver Richards equation. 
 76   
 77  kex is the chemical exchange rate constant, pA and pB are the populations of states A and B, and delta_omega is the chemical shift difference between the two states in ppm. 
 78   
 79   
 80  Links 
 81  ===== 
 82   
 83  More information on the B14 model can be found in the: 
 84   
 85      - U{relax wiki<http://wiki.nmr-relax.com/B14>}, 
 86      - U{relax manual<http://www.nmr-relax.com/manual/reduced_B14_2_site_CPMG_model.html>}, 
 87      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#B14>}. 
 88   
 89  More information on the B14 full model can be found in the: 
 90   
 91      - U{relax wiki<http://wiki.nmr-relax.com/B14_full>}, 
 92      - U{relax manual<http://www.nmr-relax.com/manual/full_B14_2_site_CPMG_model.html>}, 
 93      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#B14_full>}. 
 94   
 95   
 96  Comparison to CR72 model 
 97  ======================== 
 98   
 99  Comparison to CR72 model can be found in the: 
100   
101      - U{relax wiki<http://wiki.nmr-relax.com/CR72>}, 
102      - U{relax manual<http://www.nmr-relax.com/manual/reduced_CR72_2_site_CPMG_model.html>}, 
103      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#CR72>}. 
104   
105  Comparison to CR72 full model can be found in the: 
106   
107      - U{relax wiki<http://wiki.nmr-relax.com/CR72_full>}, 
108      - U{relax manual<http://www.nmr-relax.com/manual/full_CR72_2_site_CPMG_model.html>}, 
109      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#CR72_full>}. 
110  """ 
111   
112  # Python module imports. 
113  from numpy import arccosh, arctan2, array, cos, cosh, isfinite, log, max, power, sin, sinh, sqrt, sum 
114   
115  # Repetitive calculations (to speed up calculations). 
116  g_fact = 1/sqrt(2) 
117   
118 -def r2eff_B14(r20a=None, r20b=None, pA=None, pB=None, dw=None, kex=None, k_AB=None, k_BA=None, ncyc=None, inv_tcpmg=None, tcp=None, back_calc=None, num_points=None):
119 """Calculate the R2eff values for the CR72 model. 120 121 See the module docstring for details. 122 123 124 @keyword r20a: The R20 parameter value of state A (R2 with no exchange). 125 @type r20a: float 126 @keyword r20b: The R20 parameter value of state B (R2 with no exchange). 127 @type r20b: float 128 @keyword pA: The population of state A. 129 @type pA: float 130 @keyword pB: The population of state B. 131 @type pB: float 132 @keyword dw: The chemical exchange difference between states A and B in rad/s. 133 @type dw: float 134 @keyword kex: The kex parameter value (the exchange rate in rad/s). 135 @type kex: float 136 @keyword k_AB: The rate of exchange from site A to B (rad/s). 137 @type k_AB: float 138 @keyword k_BA: The rate of exchange from site B to A (rad/s). 139 @type k_BA: float 140 @keyword ncyc: The matrix exponential power array. The number of CPMG blocks. 141 @type ncyc: numpy int16, rank-1 array 142 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds). 143 @type inv_tcpmg: float 144 @keyword tcp: The tau_CPMG times (1 / 4.nu1). 145 @type tcp: numpy rank-1 float array 146 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies. 147 @type back_calc: numpy rank-1 float array 148 @keyword num_points: The number of points on the dispersion curve, equal to the length of the cpmg_frqs and back_calc arguments. 149 @type num_points: int 150 """ 151 152 # Catch parameter values that will result in no exchange, returning flat R2eff = R20 lines (when kex = 0.0, k_AB = 0.0). 153 if dw == 0.0 or pA == 1.0 or k_AB == 0.0: 154 back_calc[:] = array([r20a]*num_points) 155 return 156 157 # Repetitive calculations (to speed up calculations). 158 deltaR2 = r20a - r20b 159 160 # The Carver and Richards (1972) alpha_minus short notation. 161 alpha_m = deltaR2 + k_AB - k_BA 162 zeta = 2.0 * dw * alpha_m 163 Psi = alpha_m**2 + 4.0 * k_BA * k_AB - dw**2 164 165 # Repetitive calculations (to speed up calculations). 166 dw2 = dw**2 167 two_tcp = 2.0 * tcp 168 169 # Get the real and imaginary components of the exchange induced shift. 170 # Trigonometric functions faster than square roots. 171 quad_zeta2_Psi2 = (zeta**2 + Psi**2)**0.25 172 g3 = cos(0.5 * arctan2(-zeta, Psi)) * quad_zeta2_Psi2 173 g4 = sin(0.5 * arctan2(-zeta, Psi)) * quad_zeta2_Psi2 174 175 # Repetitive calculations (to speed up calculations). 176 g32 = g3**2 177 g42 = g4**2 178 179 # Time independent factors. 180 # N = oG + oE. 181 N = g3 + g4*1j 182 183 NNc = g32 + g42 184 185 # F0. 186 F0 = (dw2 + g32) / NNc 187 188 # F2. 189 F2 = (dw2 - g42) / NNc 190 191 # t1 = (-dw + g4) * (complex(-dw, -g3)) / NNc #t1. 192 193 # t2. 194 F1b = (dw + g4) * (dw - g3*1j) / NNc 195 196 # t1 + t2. 197 F1a_plus_b = (2. * dw2 + zeta*1j) / NNc 198 199 # Derived from relaxation. 200 # E0 = -2.0 * tcp * (F00R - f11R). 201 E0 = two_tcp * g3 202 203 # Catch math domain error of sinh(val > 710). 204 # This is when E0 > 710. 205 if max(E0) > 700: 206 back_calc[:] = array([r20a]*num_points) 207 return 208 209 # Derived from chemical shifts #E2 = complex(0,-2.0 * tcp * (F00I - f11I)). 210 E2 = two_tcp * g4 211 212 # Mixed term (complex) (E0 - iE2)/2. 213 E1 = (g3 - g4*1j) * tcp 214 215 # Complex. 216 v1s = F0 * sinh(E0) - F2 * sin(E2)*1j 217 218 # -2 * oG * t2. 219 v4 = F1b * (-alpha_m - g3 ) + F1b * (dw - g4)*1j 220 221 # Complex. 222 ex1c = sinh(E1) 223 224 # Off diagonal common factor. sinh fuctions. 225 v5 = (-deltaR2 + kex + dw*1j) * v1s - 2. * (v4 + k_AB * F1a_plus_b) * ex1c 226 227 # Real. The v_1c in paper. 228 v1c = F0 * cosh(E0) - F2 * cos(E2) 229 230 # Exact result for v2v3. 231 v3 = sqrt(v1c**2 - 1.) 232 233 y = power( (v1c - v3) / (v1c + v3), ncyc) 234 235 Tog = 0.5 * (1. + y) + (1. - y) * v5 / (2. * v3 * N ) 236 237 ## -1/Trel * log(LpreDyn). 238 # Rpre = (r20a + r20b + kex) / 2.0 239 240 ## Carver and Richards (1972) 241 # R2eff_CR72 = Rpre - inv_tcpmg * ncyc * arccosh(v1c.real) 242 243 ## Baldwin final. 244 # Estimate R2eff. relax_time = Trel = 1/inv_tcpmg. 245 # R2eff = R2eff_CR72 - inv_tcpmg * log(Tog.real) 246 247 # Fastest calculation. 248 R2eff = (r20a + r20b + kex) / 2.0 - inv_tcpmg * ( ncyc * arccosh(v1c.real) + log(Tog.real) ) 249 250 # Catch errors, taking a sum over array is the fastest way to check for 251 # +/- inf (infinity) and nan (not a number). 252 if not isfinite(sum(R2eff)): 253 R2eff = array([1e100]*num_points) 254 255 back_calc[:] = R2eff
256