Author: bugman
Date: Sat Jun 8 23:19:26 2013
New Revision: 19973
URL: http://svn.gna.org/viewcvs/relax?rev=19973&view=rev
Log:
Added the M61 skew model equations to the relax library.
This is the Meiboom 1961 on-resonance 2-site model for skewed populations (pA
pB).
This commit follows step 3 of the relaxation dispersion model addition
tutorial
(http://thread.gmane.org/gmane.science.nmr.relax.devel/3907).
Added:
branches/relax_disp/lib/dispersion/m61b.py
- copied, changed from r19967, branches/relax_disp/lib/dispersion/m61.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=19973&r1=19972&r2=19973&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/__init__.py (original)
+++ branches/relax_disp/lib/dispersion/__init__.py Sat Jun 8 23:19:26 2013
@@ -26,5 +26,6 @@
'cr72',
'equations',
'lm63',
- 'm61'
+ 'm61',
+ 'm61b'
]
Copied: branches/relax_disp/lib/dispersion/m61b.py (from r19967,
branches/relax_disp/lib/dispersion/m61.py)
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/m61b.py?p2=branches/relax_disp/lib/dispersion/m61b.py&p1=branches/relax_disp/lib/dispersion/m61.py&r1=19967&r2=19973&rev=19973&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/m61.py (original)
+++ branches/relax_disp/lib/dispersion/m61b.py Sat Jun 8 23:19:26 2013
@@ -21,44 +21,40 @@
###############################################################################
# Module docstring.
-"""The Meiboom (1961) 2-site fast exchange R1rho model.
+"""The Meiboom (1961) 2-site on-resonance skewed population R1rho model.
-This module is for the function, gradient and Hessian of the M61 model. The
model is named after the reference:
+This module is for the function, gradient and Hessian of the M61 skew model.
The model is named after the reference:
Meiboom S. (1961). Nuclear magnetic resonance study of the proton
transfer in water. J. Chem. Phys., 34, 375-388. (U{DOI:
10.1063/1.1700960<http://dx.doi.org/10.1063/1.1700960>}).
The equation used is:
- phi_ex * kex
- R1rho = R1rho' + sin^2(theta) * ----------------- ,
- kex^2 + omega_e^2
+ pA^2.pB.delta_omega^2.kex
+ R1rho = R1rho' + -------------------------------------- ,
+ kex^2 + pA^2.delta_omega^2 + omega_1^2
-where R1rho' is the R1rho value in the absence of exchange, theta is the
rotating frame tilt angle,
-
- phi_ex = pA * pB * delta_omega^2 ,
-
-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_e is the effective field in the rotating frame.
+where 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 = omega_e is the effective field in the rotating frame.
"""
# Python module imports.
-from math import pi, sin
+from math import pi
-def r2eff_M61(r1rho_prime=None, phi_ex=None, kex=None, theta=pi/2,
spin_lock_fields=None, back_calc=None, num_points=None):
- """Calculate the R2eff values for the M61 model.
+def r1rho_M61b(r1rho_prime=None, pA=None, dw=None, kex=None,
spin_lock_fields=None, back_calc=None, num_points=None):
+ """Calculate the R1rho values for the M61 skew model.
See the module docstring for details.
@keyword r1rho_prime: The R1rho_prime parameter value (R1rho with
no exchange).
@type r1rho_prime: 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 theta: The rotating frame tilt angle.
- @type theta: float
- @keyword spin_lock_fields: The CPMG nu1 frequencies.
+ @keyword spin_lock_fields: The spin-lock field strengths (Hz).
@type spin_lock_fields: numpy rank-1 float array
@keyword back_calc: The array for holding the back calculated
R1rho values. Each element corresponds to one of the spin-lock fields.
@type back_calc: numpy rank-1 float array
@@ -66,16 +62,23 @@
@type num_poinst: int
"""
- # Loop over the time points, back calculating the R2eff values.
+ # Loop over the dispersion points, back calculating the R1rho values.
for i in range(num_points):
# Catch zeros (to avoid pointless mathematical operations).
- if phi_ex == 0.0 or kex == 0.0:
+ if pA == 0.0 or delta_omega == 0.0 or kex == 0.0:
back_calc[i] = r1rho_prime
+ continue
+
+ # Replicated calculation.
+ pA2dw2 = pA**2 * delta_omega**2
+
+ # Denominator.
+ denom = kex**2 + pA2dw2 + (2.0*pi*spin_lock_fields[i])**2
# Avoid divide by zero.
- elif kex == 0.0 and spin_lock_fields[i] == 0.0:
+ if denom == 0.0:
back_calc[i] = 1e100
+ continue
- # The full formula.
- else:
- back_calc[i] = r1rho_prime + sin(theta)**2 * phi_ex * kex /
(kex**2 + (2.0*pi*spin_lock_fields[i])**2)
+ # R1rho calculation.
+ back_calc[i] = r1rho_prime + pA2dw2 * pB * kex / denom
r19973 - in /branches/relax_disp/lib/dispersion: __init__.py m61b.py