r20008 - /branches/relax_disp/lib/dispersion/it99.py -- June 10, 2013 - 15:22

Author: bugman
Date: Mon Jun 10 15:22:01 2013
New Revision: 20008

URL: http://svn.gna.org/viewcvs/relax?rev=20008&view=rev
Log:
Added the IT99 model equations to the relax library.

This is the Ishima and Torchia 1999 2-site model for all timescales with 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/it99.py
      - copied, changed from r20006, 
branches/relax_disp/lib/dispersion/lm63.py

Copied: branches/relax_disp/lib/dispersion/it99.py (from r20006, 
branches/relax_disp/lib/dispersion/lm63.py)
URL: 
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/it99.py?p2=branches/relax_disp/lib/dispersion/it99.py&p1=branches/relax_disp/lib/dispersion/lm63.py&r1=20006&r2=20008&rev=20008&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/lm63.py (original)
+++ branches/relax_disp/lib/dispersion/it99.py Mon Jun 10 15:22:01 2013
@@ -21,31 +21,33 @@
 
###############################################################################
 
 # Module docstring.
-"""The Luz and Meiboom (1963) 2-site fast exchange model.
+"""The Ishima and Torchia (1999) 2-site model for all time scales with 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 IT99 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. 39, 366-370 (U{DOI: 
10.1063/1.1734254<http://dx.doi.org/10.1063/1.1734254>}).
+    Ishima R. and Torchia D.A. (1999).  Estimating the time scale of 
chemical exchange of proteins from measurements of transverse relaxation 
rates in solution.  J. Biomol. NMR, 14, 369-372.  (U{DOI: 
10.1023/A:1008324025406<http://dx.doi.org/10.1023/A:1008324025406>}).
 
 The equation used is:
 
-                  phi_ex   /     4 * nu_cpmg         /     kex     \ \ 
-    R2eff = R20 + ------ * | 1 - -----------  * tanh | ----------- | | ,
-                   kex     \         kex             \ 4 * nu_cpmg / /
-
-where:
+              phi_ex * tex
+    Rex ~= ------------------- ,
+           1 + omega_a^2*tex^2
 
     phi_ex = pA * pB * delta_omega^2 ,
+    
+    omega_a^2 = sqrt(omega_1^4 + pA^2*delta_omega^4) ,
 
-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.
+    R2eff = R20 + Rex ,
+
+where tex = 1/(2kex), 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.
 """
 
 # Python module imports.
-from math import tanh
+from math import 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_IT99(r20=None, phi_ex=None, padw2=None, tex=None, cpmg_frqs=None, 
back_calc=None, num_points=None):
+    """Calculate the R2eff values for the IT99 model.
 
     See the module docstring for details.
 
@@ -54,8 +56,10 @@
     @type r20:              float
     @keyword phi_ex:        The phi_ex parameter value (pA * pB * 
delta_omega^2).
     @type phi_ex:           float
-    @keyword kex:           The kex parameter value (the exchange rate in 
rad/s).
-    @type kex:              float
+    @keyword padw2:         The pA.dw^2 parameter value.
+    @type padw2:            float
+    @keyword tex:           The tex parameter value (the time of exchange in 
s/rad).
+    @type tex:              float
     @keyword cpmg_frqs:     The CPMG nu1 frequencies.
     @type cpmg_frqs:        numpy rank-1 float array
     @keyword back_calc:     The array for holding the back calculated R2eff 
values.  Each element corresponds to one of the CPMG nu1 frequencies.
@@ -65,19 +69,28 @@
     """
 
     # Repetitive calculations (to speed up calculations).
-    rex = phi_ex / kex
-    kex_4 = 4.0 / kex
+    kex = 0.5 / tex
+    tex2 = tex**2
+    pa2dw4 = padw2**2
+
+    # The numerator.
+    numer = phi_ex * tex
 
     # Loop over the time points, back calculating the R2eff values.
     for i in range(num_points):
-        # Catch zeros.
-        if phi_ex == 0.0:
+        # Catch zeros (to avoid pointless mathematical operations).
+        if numer == 0.0:
             back_calc[i] = r20
+            continue
+
+        # Denominator.
+        omega_a2 = sqrt((2.0*pi*cpmg_frqs[i])**4 + pa2dw4
+        denom = 1.0 + omega_a2 * tex2
 
         # Avoid divide by zero.
-        elif kex == 0.0:
+        if denom == 0.0:
             back_calc[i] = 1e100
+            continue
 
-        # The full formula.
-        else:
-            back_calc[i] = r20 + rex * (1.0 - kex_4 * cpmg_frqs[i] * 
tanh(kex / (4.0 * cpmg_frqs[i])))
+        # R2eff calculation.
+        back_calc[i] = r20 + numer / demon