Author: bugman
Date: Wed Aug 28 18:43:28 2013
New Revision: 20707
URL: http://svn.gna.org/viewcvs/relax?rev=20707&view=rev
Log:
Fix for the TP02 dispersion model.
The rotating frame tilt angle for this model is calculated from the
population averaged chemical
shift and not the equal weighted average.
Modified:
branches/relax_disp/lib/dispersion/tp02.py
branches/relax_disp/target_functions/relax_disp.py
Modified: branches/relax_disp/lib/dispersion/tp02.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/tp02.py?rev=20707&r1=20706&r2=20707&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/tp02.py (original)
+++ branches/relax_disp/lib/dispersion/tp02.py Wed Aug 28 18:43:28 2013
@@ -36,7 +36,7 @@
from math import atan, cos, pi, sin
-def r1rho_TP02(r1rho_prime=None, omega=None, pA=None, pB=None, dw=None,
kex=None, theta=pi/2, R1=0.0, spin_lock_fields=None, back_calc=None,
num_points=None):
+def r1rho_TP02(r1rho_prime=None, omega=None, offset=None, pA=None, pB=None,
dw=None, kex=None, R1=0.0, spin_lock_fields=None, back_calc=None,
num_points=None):
"""Calculate the R1rho' values for the TP02 model.
See the module docstring for details. This is the Trott and Palmer
(2002) equation according to Korzhnev (J. Biomol. NMR (2003), 26, 39-48).
@@ -46,6 +46,8 @@
@type r1rho_prime: float
@keyword omega: The chemical shift for the spin in rad/s.
@type omega: float
+ @keyword offset: The spin-lock offsets for the data.
+ @type offset: numpy rank-1 float array
@keyword pA: The population of state A.
@type pA: float
@keyword pB: The population of state B.
@@ -54,8 +56,6 @@
@type dw: float
@keyword kex: The kex parameter value (the exchange rate
in rad/s).
@type kex: float
- @keyword theta: The rotating frame tilt angles. Each
element corresponds to one of the spin-lock fields.
- @type theta: numpy rank-1 float array
@keyword R1: The R1 relaxation rate.
@type R1: float
@keyword spin_lock_fields: The CPMG nu1 frequencies.
@@ -70,15 +70,30 @@
half_dw = dw / 2.0
Wa = omega - half_dw # Larmor frequency [s^-1].
Wb = omega + half_dw # Larmor frequency [s^-1].
+ kex2 = kex**2
# The numerator.
numer = pA * pB * dw**2 * kex
# Loop over the dispersion points, back calculating the R1rho values.
for i in range(num_points):
+ # We assume that A resonates at 0 [s^-1], without loss of generality.
+ Wsl = offset[i] # Larmor frequency of spin lock
[s^-1].
+ w1 = spin_lock_fields[i] * 2.0 * pi # Spin-lock field strength.
+ W = pA*Wa + pB*Wb # Pop-averaged Larmor frequency
[s^-1].
+ da = Wa - Wsl # Offset of spin-lock from A.
+ db = Wb - Wsl # Offset of spin-lock from B.
+ d = W - Wsl # Offset of spin-lock from
pop-average.
+ waeff2 = w1**2 + da**2 # Effective field at A.
+ wbeff2 = w1**2 + db**2 # Effective field at B.
+ weff2 = w1**2 + d**2 # Effective field at pop-average.
+
+ # The rotating frame flip angle.
+ theta = atan(w1 / d)
+
# Repetitive calculations (to speed up calculations).
- sin_theta2 = sin(theta[i])**2
- R1_cos_theta2 = R1 * cos(theta[i])**2
+ sin_theta2 = sin(theta)**2
+ R1_cos_theta2 = R1 * cos(theta)**2
R1rho_prime_sin_theta2 = r1rho_prime * sin_theta2
# Catch zeros (to avoid pointless mathematical operations).
@@ -86,20 +101,9 @@
back_calc[i] = R1_cos_theta2 + R1rho_prime_sin_theta2
continue
- # We assume that A resonates at 0 [s^-1], without loss of generality.
- Wsl = spin_lock_fields[i] * 2.0 * pi # Larmor frequency of spin
lock [s^-1].
- w1 = spin_lock_fields[i] * 2.0 * pi # Larmor frequency of spin
lock [s^-1].
- W = pA*Wa + pB*Wb # Pop-averaged Larmor frequency [s^-1].
- da = Wa - Wsl # Offset of spin-lock from A.
- db = Wb - Wsl # Offset of spin-lock from B.
- d = W - Wsl # Offset of spin-lock from pop-average.
- waeff2 = w1**2 + da**2 # Effective field at A.
- wbeff2 = w1**2 + db**2 # Effective field at B.
- weff2 = w1**2 + d**2 # Effective field at pop-average.
-
# Denominator.
- denom = waeff2 * wbeff2 / weff2 + kex**2
- #denom_extended = waeff2*wbeff2/weff2+kex**2-2*sin_theta2*pA*pB*dw**2
+ denom = waeff2 * wbeff2 / weff2 + kex2
+ #denom_extended = waeff2*wbeff2/weff2+kex2-2*sin_theta2*pA*pB*dw**2
# Avoid divide by zero.
if denom == 0.0:
Modified: branches/relax_disp/target_functions/relax_disp.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/target_functions/relax_disp.py?rev=20707&r1=20706&r2=20707&view=diff
==============================================================================
--- branches/relax_disp/target_functions/relax_disp.py (original)
+++ branches/relax_disp/target_functions/relax_disp.py Wed Aug 28 18:43:28
2013
@@ -965,7 +965,7 @@
dw_frq = dw[spin_index] * self.frqs[spin_index, frq_index]
# Back calculate the R1rho values.
- r1rho_TP02(r1rho_prime=R20[r20_index],
omega=self.chemical_shifts[spin_index, frq_index], pA=pA, pB=pB, dw=dw_frq,
kex=kex, theta=self.tilt_angles[spin_index, frq_index],
R1=self.r1[spin_index, frq_index], spin_lock_fields=self.spin_lock_nu1,
back_calc=self.back_calc[spin_index, frq_index],
num_points=self.num_disp_points)
+ r1rho_TP02(r1rho_prime=R20[r20_index],
omega=self.chemical_shifts[spin_index, frq_index],
offset=self.spin_lock_offsets[spin_index, frq_index], pA=pA, pB=pB,
dw=dw_frq, kex=kex, R1=self.r1[spin_index, frq_index],
spin_lock_fields=self.spin_lock_nu1, back_calc=self.back_calc[spin_index,
frq_index], num_points=self.num_disp_points)
# For all missing data points, set the back-calculated value
to the measured values so that it has no effect on the chi-squared value.
for point_index in range(self.num_disp_points):
r20707 - in /branches/relax_disp: lib/dispersion/tp02.py target_functions/relax_disp.py