Author: bugman
Date: Wed Dec 4 11:29:25 2013
New Revision: 21764
URL: http://svn.gna.org/viewcvs/relax?rev=21764&view=rev
Log:
Fixes for the 'NS MMQ 3-site' dispersion models - the evolution matrix is now
correctly constructed.
Modified:
trunk/lib/dispersion/ns_mmq_3site.py
trunk/target_functions/relax_disp.py
Modified: trunk/lib/dispersion/ns_mmq_3site.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/lib/dispersion/ns_mmq_3site.py?rev=21764&r1=21763&r2=21764&view=diff
==============================================================================
--- trunk/lib/dispersion/ns_mmq_3site.py (original)
+++ trunk/lib/dispersion/ns_mmq_3site.py Wed Dec 4 11:29:25 2013
@@ -39,7 +39,7 @@
from lib.linear_algebra.matrix_power import square_matrix_power
-def populate_matrix(matrix=None, R20A=None, R20B=None, R20C=None,
dw_AB=None, dw_BC=None, dw_AC=None, k_AB=None, k_BA=None, k_BC=None,
k_CB=None, k_AC=None, k_CA=None):
+def populate_matrix(matrix=None, R20A=None, R20B=None, R20C=None,
dw_AB=None, dw_AC=None, k_AB=None, k_BA=None, k_BC=None, k_CB=None,
k_AC=None, k_CA=None):
"""The Bloch-McConnell matrix for 3-site exchange.
@keyword matrix: The matrix to populate.
@@ -52,8 +52,6 @@
@type R20C: float
@keyword dw_AB: The combined chemical exchange difference
parameters between states A and B in rad/s. This can be any combination of
dw and dwH.
@type dw_AB: float
- @keyword dw_BC: The combined chemical exchange difference
parameters between states B and C in rad/s. This can be any combination of
dw and dwH.
- @type dw_BC: float
@keyword dw_AC: The combined chemical exchange difference
parameters between states A and C in rad/s. This can be any combination of
dw and dwH.
@type dw_AC: float
@keyword k_AB: The rate of exchange from site A to B (rad/s).
@@ -71,7 +69,7 @@
"""
# The first row.
- matrix[0, 0] = -k_AB - k_AC + 1.j*dw_BC - R20A
+ matrix[0, 0] = -k_AB - k_AC - R20A
matrix[0, 1] = k_BA
matrix[0, 2] = k_CA
@@ -86,7 +84,7 @@
matrix[2, 2] = -k_CB - k_CA + 1.j*dw_AC - R20C
-def r2eff_ns_mmq_3site_mq(M0=None, F_vector=array([1, 0, 0], float64),
m1=None, m2=None, R20A=None, R20B=None, R20C=None, pA=None, pB=None, pC=None,
dw_AB=None, dw_BC=None, dw_AC=None, dwH_AB=None, dwH_BC=None, dwH_AC=None,
k_AB=None, k_BA=None, k_BC=None, k_CB=None, k_AC=None, k_CA=None,
inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None):
+def r2eff_ns_mmq_3site_mq(M0=None, F_vector=array([1, 0, 0], float64),
m1=None, m2=None, R20A=None, R20B=None, R20C=None, pA=None, pB=None, pC=None,
dw_AB=None, dw_AC=None, dwH_AB=None, dwH_AC=None, k_AB=None, k_BA=None,
k_BC=None, k_CB=None, k_AC=None, k_CA=None, inv_tcpmg=None, tcp=None,
back_calc=None, num_points=None, power=None):
"""The 3-site numerical solution to the Bloch-McConnell equation for MQ
data.
The notation used here comes from:
@@ -122,14 +120,10 @@
@type pC: float
@keyword dw_AB: The chemical exchange difference between states
A and B in rad/s.
@type dw_AB: float
- @keyword dw_BC: The chemical exchange difference between states
B and C in rad/s.
- @type dw_BC: float
@keyword dw_AC: The chemical exchange difference between states
A and C in rad/s.
@type dw_AC: float
@keyword dwH_AB: The proton chemical exchange difference between
states A and B in rad/s.
@type dwH_AB: float
- @keyword dwH_BC: The proton chemical exchange difference between
states B and C in rad/s.
- @type dwH_BC: float
@keyword dwH_AC: The proton chemical exchange difference between
states A and C in rad/s.
@type dwH_AC: float
@keyword k_AB: The rate of exchange from site A to B (rad/s).
@@ -157,8 +151,8 @@
"""
# Populate the m1 and m2 matrices (only once per function call for
speed).
- populate_matrix(matrix=m1, R20A=R20A, R20B=R20B, R20C=R20C,
dw_AB=-dw_AB-dwH_AB, dw_BC=-dw_BC-dwH_BC, dw_AC=-dw_AC-dwH_AC, k_AB=k_AB,
k_BA=k_BA, k_BC=k_BC, k_CB=k_CB, k_AC=k_AC, k_CA=k_CA) # D+ matrix
component.
- populate_matrix(matrix=m2, R20A=R20A, R20B=R20B, R20C=R20C,
dw_AB=-dw_AB-dwH_AB, dw_BC=dw_BC-dwH_BC, dw_AC=dw_AC-dwH_AC, k_AB=k_AB,
k_BA=k_BA, k_BC=k_BC, k_CB=k_CB, k_AC=k_AC, k_CA=k_CA) # Z- matrix
component.
+ populate_matrix(matrix=m1, R20A=R20A, R20B=R20B, R20C=R20C,
dw_AB=-dw_AB-dwH_AB, dw_AC=-dw_AC-dwH_AC, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC,
k_CB=k_CB, k_AC=k_AC, k_CA=k_CA) # D+ matrix component.
+ populate_matrix(matrix=m2, R20A=R20A, R20B=R20B, R20C=R20C,
dw_AB=-dw_AB-dwH_AB, dw_AC=dw_AC-dwH_AC, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC,
k_CB=k_CB, k_AC=k_AC, k_CA=k_CA) # Z- matrix component.
# Loop over the time points, back calculating the R2eff values.
for i in range(num_points):
@@ -241,9 +235,10 @@
back_calc[i] = 1e99
else:
back_calc[i]= -inv_tcpmg * log(Mx / pA)
-
-
-def r2eff_ns_mmq_3site_sq_dq_zq(M0=None, F_vector=array([1, 0, 0], float64),
m1=None, m2=None, R20A=None, R20B=None, R20C=None, pA=None, pB=None, pC=None,
dw_AB=None, dw_BC=None, dw_AC=None, dwH_AB=None, dwH_BC=None, dwH_AC=None,
k_AB=None, k_BA=None, k_BC=None, k_CB=None, k_AC=None, k_CA=None,
inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None):
+ print back_calc
+
+
+def r2eff_ns_mmq_3site_sq_dq_zq(M0=None, F_vector=array([1, 0, 0], float64),
m1=None, m2=None, R20A=None, R20B=None, R20C=None, pA=None, pB=None, pC=None,
dw_AB=None, dw_AC=None, dwH_AB=None, dwH_AC=None, k_AB=None, k_BA=None,
k_BC=None, k_CB=None, k_AC=None, k_CA=None, inv_tcpmg=None, tcp=None,
back_calc=None, num_points=None, power=None):
"""The 3-site numerical solution to the Bloch-McConnell equation for SQ,
ZQ, and DQ data.
The notation used here comes from:
@@ -275,14 +270,10 @@
@type pC: float
@keyword dw_AB: The combined chemical exchange difference
between states A and B in rad/s. It should be set to dwH for 1H SQ data, dw
for heteronuclear SQ data, dwH-dw for ZQ data, and dwH+dw for DQ data.
@type dw_AB: float
- @keyword dw_BC: The combined chemical exchange difference
between states B and C in rad/s. It should be set to dwH for 1H SQ data, dw
for heteronuclear SQ data, dwH-dw for ZQ data, and dwH+dw for DQ data.
- @type dw_BC: float
@keyword dw_AC: The combined chemical exchange difference
between states A and C in rad/s. It should be set to dwH for 1H SQ data, dw
for heteronuclear SQ data, dwH-dw for ZQ data, and dwH+dw for DQ data.
@type dw_AB: float
@keyword dwH_AB: Unused - this is simply to match the
r2eff_mmq_3site_mq() function arguments.
@type dwH_AC: float
- @keyword dwH_BC: Unused - this is simply to match the
r2eff_mmq_3site_mq() function arguments.
- @type dwH_BC: float
@keyword dwH_AC: Unused - this is simply to match the
r2eff_mmq_3site_mq() function arguments.
@type dwH_AB: float
@keyword k_AB: The rate of exchange from site A to B (rad/s).
@@ -302,8 +293,8 @@
"""
# Populate the m1 and m2 matrices (only once per function call for
speed).
- populate_matrix(matrix=m1, R20A=R20A, R20B=R20B, R20C=R20C, dw_AB=dw_AB,
dw_BC=dw_BC, dw_AC=dw_AC, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC, k_CB=k_CB,
k_AC=k_AC, k_CA=k_CA)
- populate_matrix(matrix=m2, R20A=R20A, R20B=R20B, R20C=R20C,
dw_AB=-dw_AB, dw_BC=-dw_BC, dw_AC=-dw_AC, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC,
k_CB=k_CB, k_AC=k_AC, k_CA=k_CA)
+ populate_matrix(matrix=m1, R20A=R20A, R20B=R20B, R20C=R20C, dw_AB=dw_AB,
dw_AC=dw_AC, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC, k_CB=k_CB, k_AC=k_AC, k_CA=k_CA)
+ populate_matrix(matrix=m2, R20A=R20A, R20B=R20B, R20C=R20C,
dw_AB=-dw_AB, dw_AC=-dw_AC, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC, k_CB=k_CB,
k_AC=k_AC, k_CA=k_CA)
# Loop over the time points, back calculating the R2eff values.
for i in range(num_points):
Modified: trunk/target_functions/relax_disp.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/target_functions/relax_disp.py?rev=21764&r1=21763&r2=21764&view=diff
==============================================================================
--- trunk/target_functions/relax_disp.py (original)
+++ trunk/target_functions/relax_disp.py Wed Dec 4 11:29:25 2013
@@ -495,7 +495,7 @@
return chi2_sum
- def calc_ns_mmq_3site_chi2(self, R20A=None, R20B=None, R20C=None,
dw_AB=None, dw_BC=None, dw_AC=None, dwH_AB=None, dwH_BC=None, dwH_AC=None,
pA=None, pB=None, kex_AB=None, kex_BC=None, kex_AC=None):
+ def calc_ns_mmq_3site_chi2(self, R20A=None, R20B=None, R20C=None,
dw_AB=None, dw_BC=None, dwH_AB=None, dwH_BC=None, pA=None, pB=None,
kex_AB=None, kex_BC=None, kex_AC=None):
"""Calculate the chi-squared value for the 'NS MMQ 3-site' models.
@keyword R20A: The R2 value for state A in the absence of
exchange.
@@ -526,12 +526,15 @@
# Once off parameter conversions.
pC = 1.0 - pA - pB
- k_BA = pA * kex_AB
- k_AB = pB * kex_AB
- k_CA = pA * kex_AC
- k_AC = pC * kex_AC
- k_CB = pB * kex_BC
- k_BC = pC * kex_BC
+ pA_pB = pA + pB
+ pA_pC = pA + pC
+ pB_pC = pB + pC
+ k_BA = pA * kex_AB / pA_pB
+ k_AB = pB * kex_AB / pA_pB
+ k_CB = pB * kex_BC / pB_pC
+ k_BC = pC * kex_BC / pB_pC
+ k_CA = pA * kex_AC / pA_pC
+ k_AC = pC * kex_AC / pA_pC
dw_AC = dw_AB + dw_BC
dwH_AC = dwH_AB + dwH_BC
@@ -554,49 +557,38 @@
# Convert dw from ppm to rad/s.
dw_AB_frq = dw_AB[spin_index] *
self.frqs[exp_index][spin_index][frq_index]
- dw_BC_frq = dw_BC[spin_index] *
self.frqs[exp_index][spin_index][frq_index]
dw_AC_frq = dw_AC[spin_index] *
self.frqs[exp_index][spin_index][frq_index]
dwH_AB_frq = dwH_AB[spin_index] *
self.frqs_H[exp_index][spin_index][frq_index]
- dwH_BC_frq = dwH_BC[spin_index] *
self.frqs_H[exp_index][spin_index][frq_index]
dwH_AC_frq = dwH_AC[spin_index] *
self.frqs_H[exp_index][spin_index][frq_index]
# Alias the dw frequency combinations.
aliased_dwH_AB = 0.0
- aliased_dwH_BC = 0.0
aliased_dwH_AC = 0.0
if self.exp_types[exp_index] == EXP_TYPE_CPMG_SQ:
aliased_dw_AB = dw_AB_frq
- aliased_dw_BC = dw_BC_frq
aliased_dw_AC = dw_AC_frq
elif self.exp_types[exp_index] ==
EXP_TYPE_CPMG_PROTON_SQ:
aliased_dw_AB = dwH_AB_frq
- aliased_dw_BC = dwH_BC_frq
aliased_dw_AC = dwH_AC_frq
elif self.exp_types[exp_index] == EXP_TYPE_CPMG_DQ:
aliased_dw_AB = dw_AB_frq + dwH_AB_frq
- aliased_dw_BC = dw_BC_frq + dwH_BC_frq
aliased_dw_AC = dw_AC_frq + dwH_AC_frq
elif self.exp_types[exp_index] == EXP_TYPE_CPMG_ZQ:
aliased_dw_AB = dw_AB_frq - dwH_AB_frq
- aliased_dw_BC = dw_BC_frq - dwH_BC_frq
aliased_dw_AC = dw_AC_frq - dwH_AC_frq
elif self.exp_types[exp_index] == EXP_TYPE_CPMG_MQ:
aliased_dw_AB = dw_AB_frq
- aliased_dw_BC = dw_BC_frq
aliased_dw_AC = dw_AC_frq
aliased_dwH_AB = dwH_AB_frq
- aliased_dwH_BC = dwH_BC_frq
aliased_dwH_AC = dwH_AC_frq
elif self.exp_types[exp_index] ==
EXP_TYPE_CPMG_PROTON_MQ:
aliased_dw_AB = dwH_AB_frq
- aliased_dw_BC = dwH_BC_frq
aliased_dw_AC = dwH_AC_frq
aliased_dwH_AB = dw_AB_frq
- aliased_dwH_BC = dw_BC_frq
aliased_dwH_AC = dw_AC_frq
# Back calculate the R2eff values for each experiment
type.
- self.r2eff_mmq[exp_index](M0=self.M0, m1=self.m1,
m2=self.m2, R20A=R20A[r20_index], R20B=R20B[r20_index], R20C=R20C[r20_index],
pA=pA, pB=pB, pC=pC, dw_AB=aliased_dw_AB, dw_BC=aliased_dw_BC,
dw_AC=aliased_dw_AC, dwH_AB=aliased_dwH_AB, dwH_BC=aliased_dwH_BC,
dwH_AC=aliased_dwH_AC, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC, k_CB=k_CB, k_AC=k_AC,
k_CA=k_CA, inv_tcpmg=self.inv_relax_times[exp_index][frq_index],
tcp=self.tau_cpmg[exp_index][frq_index],
back_calc=self.back_calc[exp_index][spin_index][frq_index],
num_points=self.num_disp_points[exp_index][frq_index],
power=self.power[exp_index][frq_index])
+ self.r2eff_mmq[exp_index](M0=self.M0, m1=self.m1,
m2=self.m2, R20A=R20A[r20_index], R20B=R20B[r20_index], R20C=R20C[r20_index],
pA=pA, pB=pB, pC=pC, dw_AB=aliased_dw_AB, dw_AC=aliased_dw_AC,
dwH_AB=aliased_dwH_AB, dwH_AC=aliased_dwH_AC, k_AB=k_AB, k_BA=k_BA,
k_BC=k_BC, k_CB=k_CB, k_AC=k_AC, k_CA=k_CA,
inv_tcpmg=self.inv_relax_times[exp_index][frq_index],
tcp=self.tau_cpmg[exp_index][frq_index],
back_calc=self.back_calc[exp_index][spin_index][frq_index],
num_points=self.num_disp_points[exp_index][frq_index],
power=self.power[exp_index][frq_index])
# 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[exp_index][frq_index]):
r21764 - in /trunk: lib/dispersion/ns_mmq_3site.py target_functions/relax_disp.py