lib.dispersion.cr72

Source Code for Module lib.dispersion.cr72

1 ############################################################################### 2 # # 3 # Copyright (C) 2009 Sebastien Morin # 4 # Copyright (C) 2013-2014 Edward d'Auvergne # 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 Carver and Richards (1972) 2-site all time scale exchange U{CR72<http://wiki.nmr-relax.com/CR72>} and U{CR72 full<http://wiki.nmr-relax.com/CR72_full>} models. 25 26 Description 27 =========== 28 29 This module is for the function, gradient and Hessian of the U{CR72<http://wiki.nmr-relax.com/CR72>} and U{CR72 full<http://wiki.nmr-relax.com/CR72_full>} models. 30 31 32 References 33 ========== 34 35 The model is named after the reference: 36 37 - Carver, J. P. and Richards, R. E. (1972). General 2-site solution for chemical exchange produced dependence of T2 upon Carr-Purcell pulse separation. I{J. Magn. Reson.}, B{6}, 89-105. (U{DOI: 10.1016/0022-2364(72)90090-X<http://dx.doi.org/10.1016/0022-2364(72)90090-X>}). 38 39 40 Equations 41 ========= 42 43 The equation used is:: 44 45 R2eff = 1/2 [ R2A0 + R2B0 + kex - 2.nu_cpmg.cosh^-1 (D+.cosh(eta+) - D-.cos(eta-)) ] , 46 47 where:: 48 49 1 / Psi + 2delta_omega^2 \ 50 D+/- = - | +/-1 + -------------------- | , 51 2 \ sqrt(Psi^2 + zeta^2) / 52 53 1 54 eta+/- = 2^(-3/2) . -------- sqrt(+/-Psi + sqrt(Psi^2 + zeta^2)) , 55 nu_cpmg 56 57 Psi = (R2A0 - R2B0 - pA.kex + pB.kex)^2 - delta_omega^2 + 4pA.pB.kex^2 , 58 59 zeta = 2delta_omega (R2A0 - R2B0 - pA.kex + pB.kex). 60 61 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. 62 63 64 CR72 model 65 ---------- 66 67 Importantly for the implementation of this model, it is assumed that R2A0 and R2B0 are identical. This simplifies some of the equations to:: 68 69 R2eff = R20 + kex/2 - nu_cpmg.cosh^-1 (D+.cosh(eta+) - D-.cos(eta-) , 70 71 where:: 72 73 Psi = kex^2 - delta_omega^2 , 74 75 zeta = -2delta_omega (pA.kex - pB.kex). 76 77 78 Links 79 ===== 80 81 More information on the CR72 model can be found in the: 82 83 - U{relax wiki<http://wiki.nmr-relax.com/CR72>}, 84 - U{relax manual<http://www.nmr-relax.com/manual/reduced_CR72_2_site_CPMG_model.html>}, 85 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#CR72>}. 86 87 More information on the CR72 full model can be found in the: 88 89 - U{relax wiki<http://wiki.nmr-relax.com/CR72_full>}, 90 - U{relax manual<http://www.nmr-relax.com/manual/full_CR72_2_site_CPMG_model.html>}, 91 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#CR72_full>}. 92 """ 93 94 # Python module imports. 95 from numpy import arccosh, array, cos, cosh, isfinite, min, max, sqrt, sum 96 97 # Repetitive calculations (to speed up calculations). 98 eta_scale = 2.0**(-3.0/2.0) 99

100 -def r2eff_CR72(r20a=None, r20b=None, pA=None, dw=None, kex=None, cpmg_frqs=None, back_calc=None, num_points=None):

101 """Calculate the R2eff values for the CR72 model. 102 103 See the module docstring for details. 104 105 106 @keyword r20a: The R20 parameter value of state A (R2 with no exchange). 107 @type r20a: float 108 @keyword r20b: The R20 parameter value of state B (R2 with no exchange). 109 @type r20b: float 110 @keyword pA: The population of state A. 111 @type pA: float 112 @keyword dw: The chemical exchange difference between states A and B in rad/s. 113 @type dw: float 114 @keyword kex: The kex parameter value (the exchange rate in rad/s). 115 @type kex: float 116 @keyword cpmg_frqs: The CPMG nu1 frequencies. 117 @type cpmg_frqs: numpy rank-1 float array 118 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies. 119 @type back_calc: numpy rank-1 float array 120 @keyword num_points: The number of points on the dispersion curve, equal to the length of the cpmg_frqs and back_calc arguments. 121 @type num_points: int 122 """ 123 124 # Catch parameter values that will result in no exchange, returning flat R2eff = R20 lines (when kex = 0.0, k_AB = 0.0). 125 if dw == 0.0 or pA == 1.0 or kex == 0.0: 126 back_calc[:] = array([r20a]*num_points) 127 return 128 129 # The B population. 130 pB = 1.0 - pA 131 132 # Repetitive calculations (to speed up calculations). 133 dw2 = dw**2 134 r20_kex = (r20a + r20b + kex) / 2.0 135 k_BA = pA * kex 136 k_AB = pB * kex 137 138 # The Psi and zeta values. 139 if r20a != r20b: 140 fact = r20a - r20b - k_BA + k_AB 141 Psi = fact**2 - dw2 + 4.0*pA*pB*kex**2 142 zeta = 2.0*dw * fact 143 else: 144 Psi = kex**2 - dw2 145 zeta = -2.0*dw * (k_BA - k_AB) 146 147 # More repetitive calculations. 148 sqrt_psi2_zeta2 = sqrt(Psi**2 + zeta**2) 149 150 # The D+/- values. 151 D_part = (Psi + 2.0*dw2) / sqrt_psi2_zeta2 152 Dpos = 0.5 * (1.0 + D_part) 153 Dneg = 0.5 * (-1.0 + D_part) 154 155 # Partial eta+/- values. 156 etapos = eta_scale * sqrt(Psi + sqrt_psi2_zeta2) / cpmg_frqs 157 etaneg = eta_scale * sqrt(-Psi + sqrt_psi2_zeta2) / cpmg_frqs 158 159 # Catch math domain error of cosh(val > 710). 160 # This is when etapos > 710. 161 if max(etapos) > 700: 162 back_calc[:] = array([r20a]*num_points) 163 return 164 165 # The arccosh argument - catch invalid values. 166 fact = Dpos * cosh(etapos) - Dneg * cos(etaneg) 167 if min(fact) < 1.0: 168 back_calc[:] = array([r20_kex]*num_points) 169 return 170 171 # Calculate R2eff. 172 R2eff = r20_kex - cpmg_frqs * arccosh( fact ) 173 174 # Catch errors, taking a sum over array is the fastest way to check for 175 # +/- inf (infinity) and nan (not a number). 176 if not isfinite(sum(R2eff)): 177 R2eff = array([1e100]*num_points) 178 179 back_calc[:] = R2eff

180