Author: bugman
Date: Fri May 31 10:39:48 2013
New Revision: 19814
URL: http://svn.gna.org/viewcvs/relax?rev=19814&view=rev
Log:
Added the CR72 model equations to the relax library.
This is for the Carver and Richards 1972 2-site exchange model covering all
time scales.
Added:
branches/relax_disp/lib/dispersion/cr72.py
- copied, changed from r19805,
branches/relax_disp/lib/dispersion/lm63.py
Modified:
branches/relax_disp/lib/dispersion/__init__.py
Modified: branches/relax_disp/lib/dispersion/__init__.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/__init__.py?rev=19814&r1=19813&r2=19814&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/__init__.py (original)
+++ branches/relax_disp/lib/dispersion/__init__.py Fri May 31 10:39:48 2013
@@ -23,5 +23,7 @@
"""The relax-lib NMR package - a library of functions for relaxation
dispersion."""
__all__ = [
- 'equations'
+ 'cr72',
+ 'equations',
+ 'lm63'
]
Copied: branches/relax_disp/lib/dispersion/cr72.py (from r19805,
branches/relax_disp/lib/dispersion/lm63.py)
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/cr72.py?p2=branches/relax_disp/lib/dispersion/cr72.py&p1=branches/relax_disp/lib/dispersion/lm63.py&r1=19805&r2=19814&rev=19814&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/lm63.py (original)
+++ branches/relax_disp/lib/dispersion/cr72.py Fri May 31 10:39:48 2013
@@ -21,39 +21,56 @@
###############################################################################
# Module docstring.
-"""The Luz and Meiboom (1963) 2-site fast exchange model.
+"""The Carver and Richards (1972) 2-site all time scale exchange model (pA >
pB).
-This module is for the function, gradient and Hessian of the LM63 model.
The model is named after the reference:
+This module is for the function, gradient and Hessian of the CR72 model.
The model is named after the reference:
- 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, J. Chem. Phys. 1963, 39, 366-370 (U{DOI:
10.1063/1.1734254<http://dx.doi.org/10.1063/1.1734254>}).
+ Carver, J. P. and Richards, R. E. (1972). General 2-site solution for
chemical exchange produced dependence of T2 upon Carr-Purcell pulse
separation. J. Magn. Reson., 6, 89-105. (U{DOI:
10.1016/0022-2364(72)90090-X<http://dx.doi.org/10.1016/0022-2364(72)90090-X>}).
The equation used is:
- phi_ex / 4 * nu_cpmg / kex \ \
- R2eff = R20 + ------ * | 1 - ----------- * tanh | ----------- | | ,
- kex \ kex \ 4 * nu_cpmg / /
+ R2eff = 1/2 [ R2A0 + R2B0 + kex - 2.nu_cpmg.cosh^-1 (D+.cosh(eta+) -
D-.cos(eta-) ] ,
+
+where:
+ 1 / Psi + 2delta_omega^2 \
+ D+/- = - | +/-1 + -------------------- | ,
+ 2 \ sqrt(Psi^2 + zeta^2) /
+
+ 2^(2/3)
+ eta+/- = ------- sqrt(+/-Psi + sqrt(Psi^2 + zeta^2)) ,
+ nu_cpmg
+
+ Psi = (R2A0 - R2B0 - pA.kex + pB.kex)^2 - delta_omega^2 + 4pA.pB.kex^2 ,
+
+ zeta = 2delta_omega (R2A0 - R2B0 - pA.kex + pB.kex).
+
+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. Importantly for the implementation of this model, it is
assumed that R2A0 and R2B0 are identical. This simplifies some of the
equations to:
+
+ R2eff = R20 + kex/2 - nu_cpmg.cosh^-1 (D+.cosh(eta+) - D-.cos(eta-) ,
where:
- phi_ex = pA * pB * delta_omega^2 ,
+ Psi = kex^2 - delta_omega^2 ,
-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.
+ zeta = -2delta_omega (pA.kex - pB.kex).
"""
# Python module imports.
-from math import tanh
+from math import arccosh, cos, cosh, sqrt
-def r2eff_LM63(r20=None, phi_ex=None, kex=None, cpmg_frqs=None,
back_calc=None, num_points=None):
- """Calculate the R2eff values for the LM63 model.
+def r2eff_CR72(r20=None, pA=None, dw=None, kex=None, cpmg_frqs=None,
back_calc=None, num_points=None):
+ """Calculate the R2eff values for the CR72 model.
See the module docstring for details.
- @keyword r20: The R20 parameter value (R2 with no exchange).
+ @keyword r20: The R20 parameter value (R2 with no exchange).
It is assumed that R2A0 and R2B0 are identical.
@type r20: float
- @keyword phi_ex: The phi_ex parameter value (pA * pB *
delta_omega^2).
- @type phi_ex: float
+ @keyword pA: The population of state A.
+ @type pA: float
+ @keyword dw: The chemical exchange difference between states
A and B in rad/s.
+ @type dw: float
@keyword kex: The kex parameter value (the exchange rate in
rad/s).
@type kex: float
@keyword cpmg_frqs: The CPMG nu1 frequencies.
@@ -64,16 +81,36 @@
@type num_poinst: int
"""
+ # The B population.
+ pB = 1.0 - pA
+
+ # Repetitive calculations (to speed up calculations).
+ dw2 = dw**2
+
+ # The Psi value.
+ Psi = kex**2 - dw2
+
+ # The zeta value.
+ zeta = -2.0*dw * (pA*kex - pB*kex)
+
+ # More repetitive calculations.
+ sqrt_psi2_zeta2 = sqrt(Psi**2 + zeta**2)
+
+ # The D+/- values.
+ D_part = (Psi + 2.0*dw2) / sqrt_psi2_zeta2
+ Dpos = 0.5 * (1.0 + D_part)
+ Dneg = 0.5 * (-1.0 + D_part)
+
+ # Partial eta+/- values.
+ eta_scale = 2.0**(2.0/3.0)
+ etapos_part = eta_scale * sqrt(Psi + sqrt_psi2_zeta2)
+ etaneg_part = eta_scale * sqrt(-Psi + sqrt_psi2_zeta2)
+
# Loop over the time points, back calculating the R2eff values.
for i in range(num_points):
- # Catch zeros.
- if phi_ex == 0.0:
- back_calc[i] = r20
-
- # Avoid divide by zero.
- elif kex == 0.0:
- back_calc[i] = 1e100
+ # The full eta+/- values.
+ etapos = etapos_part / cpmg_frqs[i]
+ etaneg = etaneg_part / cpmg_frqs[i]
# The full formula.
- else:
- back_calc[i] = r20 + phi_ex / kex * (1.0 - (4.0 * cpmg_frqs[i] /
kex) * tanh(kex / (4 * cpmg_frqs[i])))
+ back_calc[i] = r20 + 0.5*kex - cpmg_frqs[i] * arccosh(Dpos *
cosh(etapos) - Dneg * cos(etaneg))
r19814 - in /branches/relax_disp/lib/dispersion: __init__.py cr72.py