Author: bugman
Date: Tue Oct 15 19:45:45 2013
New Revision: 21127
URL: http://svn.gna.org/viewcvs/relax?rev=21127&view=rev
Log:
Added the 'MQ CR72' R2eff calculating function to the relax library.
This is the Carver and Richards (1972) 2-site model expanded for MQ CPMG data
by Korzhnev et al.,
2004.
This follows the tutorial for adding relaxation dispersion models at:
http://wiki.nmr-relax.com/Tutorial_for_adding_relaxation_dispersion_models_to_relax#The_relax_library.
The corresponding target function was updated to input the correct arguments.
Added:
branches/relax_disp/lib/dispersion/mq_cr72.py
- copied, changed from r21124,
branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py
Modified:
branches/relax_disp/lib/dispersion/__init__.py
branches/relax_disp/target_functions/relax_disp.py
Modified: branches/relax_disp/lib/dispersion/__init__.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/__init__.py?rev=21127&r1=21126&r2=21127&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/__init__.py (original)
+++ branches/relax_disp/lib/dispersion/__init__.py Tue Oct 15 19:45:45 2013
@@ -30,6 +30,7 @@
'lm63_3site',
'm61',
'm61b',
+ 'mq_cr72',
'mq_ns_cpmg_2site',
'ns_cpmg_2site_3d',
'ns_cpmg_2site_expanded',
Copied: branches/relax_disp/lib/dispersion/mq_cr72.py (from r21124,
branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py)
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/mq_cr72.py?p2=branches/relax_disp/lib/dispersion/mq_cr72.py&p1=branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py&r1=21124&r2=21127&rev=21127&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py (original)
+++ branches/relax_disp/lib/dispersion/mq_cr72.py Tue Oct 15 19:45:45 2013
@@ -1,7 +1,5 @@
###############################################################################
#
#
-# Copyright (C) 2013 Mathilde Lescanne
#
-# Copyright (C) 2013 Dominique Marion
#
# Copyright (C) 2013 Edward d'Auvergne
#
#
#
# This file is part of the program relax (http://www.nmr-relax.com).
#
@@ -22,69 +20,24 @@
###############################################################################
# Module docstring.
-"""This function performs a numerical fit of 2-site Bloch-McConnell
equations for MQ CPMG-type experiments.
+"""This function performs a numerical fit of CR72 model extended for MQ CPMG
data.
-The function uses an explicit matrix that contains relaxation, exchange and
chemical shift terms. It does the 180deg pulses in the CPMG train. The
approach of Bloch-McConnell can be found in chapter 3.1 of Palmer, A. G. Chem
Rev 2004, 104, 3623-3640.
+The Carver and Richards (1972) 2-site model for all times scales was
extended for multiple quantum (MQ) CPMG data by:
-This is the model of the numerical solution for the 2-site Bloch-McConnell
equations for multi-quantum CPMG-type data. It originates as the m1 and m2
matrices and the fp0() optimization function from the fitting_main_kex.py
script from Mathilde Lescanne and Dominique Marion
(https://gna.org/task/?7712#comment7 and the files attached in that comment).
+ - Korzhnev, D. M., Kloiber, K., Kanelis, V., Tugarinov, V., and Kay, L.
E. (2004). Probing slow dynamics in high molecular weight proteins by
methyl-TROSY NMR spectroscopy: Application to a 723-residue enzyme. J. Am.
Chem. Soc., 126, 3964-3973. (U{DOI:
10.1021/ja039587i<http://dx.doi.org/10.1021/ja039587i>}).
"""
-# Dependency check module.
-import dep_check
-
# Python module imports.
-from math import fabs, log
-from numpy import array, conj, dot, float64
-from numpy.linalg import matrix_power
-
-# relax module imports.
-from lib.dispersion.ns_matrices import rcpmg_3d
-from lib.float import isNaN
-from lib.linear_algebra.matrix_exponential import matrix_exponential
+from math import acosh, cos, cosh, log, sin
+from numpy import sqrt
-def populate_matrix(matrix=None, r20=None, dw=None, dwH=None, k_AB=None,
k_BA=None):
- """Matrix for HMQC experiments.
-
- This corresponds to matrix m1 and m2 matrices of equation 2.2 from:
-
- - Korzhnev, D. M., Kloiber, K., Kanelis, V., Tugarinov, V., and Kay,
L. E. (2004). Probing slow dynamics in high molecular weight proteins by
methyl-TROSY NMR spectroscopy: Application to a 723-residue enzyme. J. Am.
Chem. Soc., 126, 3964-3973. (U{DOI:
10.1021/ja039587i<http://dx.doi.org/10.1021/ja039587i>}).
-
- @keyword matrix: The matrix to populate.
- @type matrix: numpy rank-2, 2D complex64 array
- @keyword r20: The R2 value in the absence of exchange.
- @type r20: float
- @keyword dw: The chemical exchange difference between states
A and B in rad/s.
- @type dw: float
- @keyword dwH: The proton chemical exchange difference between
states A and B in rad/s.
- @type dwH: float
- @keyword k_AB: The rate of exchange from site A to B (rad/s).
- @type k_AB: float
- @keyword k_BA: The rate of exchange from site B to A (rad/s).
- @type k_BA: float
- """
-
- # Fill in the elements.
- matrix[0, 0] = -k_AB - r20
- matrix[0, 1] = k_BA
- matrix[1, 0] = k_AB
- matrix[1, 1] = -k_BA - 1.j*(dwH + dw) - r20
-
-
-def r2eff_mq_ns_cpmg_2site(M0=None, F_vector=array([1, 0], float64),
m1=None, m2=None, r20=None, pA=None, pB=None, dw=None, dwH=None, k_AB=None,
k_BA=None, inv_tcpmg=None, tcp=None, back_calc=None, num_points=None,
power=None):
- """The 2-site numerical solution to the Bloch-McConnell equation.
+def r2eff_mq_cr72(r20=None, pA=None, pB=None, dw=None, dwH=None, kex=None,
k_AB=None, k_BA=None, cpmg_frqs=None, tcp=None, back_calc=None,
num_points=None, power=None):
+ """The CR72 model extended to MQ CPMG data.
This function calculates and stores the R2eff values.
- @keyword M0: This is a vector that contains the initial
magnetizations corresponding to the A and B state transverse magnetizations.
- @type M0: numpy float64, rank-1, 7D array
- @keyword F_vector: The observable magnitisation vector. This
defaults to [1, 0] for X observable magnitisation.
- @type F_vector: numpy rank-1, 2D float64 array
- @keyword m1: A complex numpy matrix to be populated.
- @type m1: numpy rank-2, 2D complex64 array
- @keyword m2: A complex numpy matrix to be populated.
- @type m2: numpy rank-2, 2D complex64 array
@keyword r20: The R2 value in the absence of exchange.
@type r20: float
@keyword pA: The population of state A.
@@ -95,12 +48,14 @@
@type dw: float
@keyword dwH: The proton chemical exchange difference between
states A and B in rad/s.
@type dwH: float
+ @keyword kex: The kex parameter value (the exchange rate in
rad/s).
+ @type kex: float
@keyword k_AB: The rate of exchange from site A to B (rad/s).
@type k_AB: float
@keyword k_BA: The rate of exchange from site B to A (rad/s).
@type k_BA: float
- @keyword inv_tcpmg: The inverse of the total duration of the CPMG
element (in inverse seconds).
- @type inv_tcpmg: float
+ @keyword cpmg_frqs: The CPMG nu1 frequencies.
+ @type cpmg_frqs: numpy rank-1 float array
@keyword tcp: The tau_CPMG times (1 / 4.nu1).
@type tcp: 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.
@@ -111,75 +66,69 @@
@type power: numpy int16, rank-1 array
"""
- # Populate the m1 and m2 matrices (only once per function call for
speed).
- populate_matrix(matrix=m1, r20=r20, dw=dw, dwH=dwH, k_AB=k_AB, k_BA=k_BA)
- populate_matrix(matrix=m2, r20=r20, dw=-dw, dwH=dwH, k_AB=k_AB,
k_BA=k_BA)
+ # Repetitive calculations (to speed up calculations).
+ dw2 = dw**2
+ r20_kex = r20 + kex/2.0
+ pApBkex2 = k_AB * k_BA
+ sqrt_pApBkex2 = sqrt(pApBkex2)
+ isqrt_pApBkex2 = 1.j*sqrt_pApBkex2
+ sqrt_pApB = sqrt(pA*pB)
+
+ # The d+/- values.
+ d = dw + dwH
+ dpos = d + 1.j*kex
+ dneg = d - 1.j*kex
+
+ # The z+/- values.
+ z = dw - dwH
+ zpos = z + 1.j*kex
+ zneg = z - 1.j*kex
+
+ # The Psi and zeta values.
+ fact = 1.j*dwH + k_BA - k_AB
+ Psi = fact**2 - dw2 + 4.0*pApBkex2
+ zeta = -2.0*dw * fact
+
+ # 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 = sqrt(2.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):
- # The M1 and M2 matrices.
- M1 = matrix_exponential(m1*tcp[i])
- M2 = matrix_exponential(m2*tcp[i])
+ # Alias delta.
+ delta = tcp[i]
+ # The full eta+/- values.
+ etapos = etapos_part * delta
+ etaneg = etaneg_part * delta
- # The complex conjugates M1* and M2*
- M1_star = conj(M1)
- M2_star = conj(M2)
+ # The mD value.
+ mD = isqrt_pApBkex2 / (dpos * zpos) * (zpos +
2.0*dw*sin(zpos*delta)/sin((dpos + zpos)*delta))
- # Repetitive dot products (minimised for speed).
- M1_M2 = dot(M1, M2)
- M2_M1 = dot(M2, M1)
- M1_M2_M2_M1 = dot(M1_M2, M2_M1)
- M2_M1_M1_M2 = dot(M2_M1, M1_M2)
- M1_M2_star = dot(M1_star, M2_star)
- M2_M1_star = dot(M2_star, M1_star)
- M1_M2_M2_M1_star = dot(M1_M2_star, M2_M1_star)
- M2_M1_M1_M2_star = dot(M2_M1_star, M1_M2_star)
+ # The mZ value.
+ mZ = -isqrt_pApBkex2 / (dneg * zneg) * (dneg -
2.0*dw*sin(dneg*delta)/sin((dneg + zneg)*delta))
- # Matrices for even n.
- if power[i] % 2 == 0:
- # The power factor (only calculate once).
- fact = power[i] / 2
+ # The Q value.
+ Q = 1 - mD**2 + mD*mZ - mZ**2 + 0.5*(mD + mZ)*sqrt_pApB
+ Q = Q.real
- # (M1.M2.M2.M1)^(n/2)
- A = matrix_power(M1_M2_M2_M1, fact)
+ # Part of the equation (catch values < 1 to prevent math domain
errors).
+ part = Dpos * cosh(etapos) - Dneg * cos(etaneg)
+ if part < 1.0:
+ back_calc[i] = 1e100
+ continue
- # (M2*.M1*.M1*.M2*)^(n/2)
- B = matrix_power(M2_M1_M1_M2_star, fact)
+ # The first eigenvalue.
+ lambda1 = r20_kex - cpmg_frqs[i] * acosh(part)
- # (M2.M1.M1.M2)^(n/2)
- C = matrix_power(M2_M1_M1_M2, fact)
-
- # (M1*.M2*.M2*.M1*)^(n/2)
- D = matrix_power(M1_M2_M2_M1_star, fact)
-
- # Matrices for odd n.
- else:
- # The power factor (only calculate once).
- fact = (power[i] - 1) / 2
-
- # (M1.M2.M2.M1)^((n-1)/2).M1.M2
- A = matrix_power(M1_M2_M2_M1, fact)
- A = dot(A, M1_M2)
-
- # (M1*.M2*.M2*.M1*)^((n-1)/2).M1*.M2*
- B = matrix_power(M1_M2_M2_M1_star, fact)
- B = dot(B, M1_M2_star)
-
- # (M2.M1.M1.M2)^((n-1)/2).M2.M1
- C = matrix_power(M2_M1_M1_M2, fact)
- C = dot(C, M2_M1)
-
- # (M2*.M1*.M1*.M2*)^((n-1)/2).M2*.M1*
- D = matrix_power(M2_M1_M1_M2_star, fact)
- D = dot(D, M2_M1_star)
-
- # The next lines calculate the R2eff using a two-point
approximation, i.e. assuming that the decay is mono-exponential.
- A_B = dot(A, B)
- C_D = dot(C, D)
- Mx = dot(dot(F_vector, (A_B + C_D)), M0)
- Mx = Mx.real / (2.0 * pA)
- if Mx <= 0.0 or isNaN(Mx):
- back_calc[i] = 1e99
- else:
- back_calc[i]= -inv_tcpmg * log(Mx)
+ # The full formula.
+ back_calc[i] = lambda1.real - cpmg_frqs[i] * log(Q) * power[i]
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=21127&r1=21126&r2=21127&view=diff
==============================================================================
--- branches/relax_disp/target_functions/relax_disp.py (original)
+++ branches/relax_disp/target_functions/relax_disp.py Tue Oct 15 19:45:45
2013
@@ -224,7 +224,7 @@
self.spin_lock_omega1_squared = self.spin_lock_omega1 ** 2
# The inverted relaxation delay.
- if model in [MODEL_MQ_CR72, MODEL_MQ_NS_CPMG_2SITE,
MODEL_NS_CPMG_2SITE_3D, MODEL_NS_CPMG_2SITE_3D_FULL,
MODEL_NS_CPMG_2SITE_EXPANDED, MODEL_NS_CPMG_2SITE_STAR,
MODEL_NS_CPMG_2SITE_STAR_FULL, MODEL_NS_R1RHO_2SITE]:
+ if model in [MODEL_MQ_NS_CPMG_2SITE, MODEL_NS_CPMG_2SITE_3D,
MODEL_NS_CPMG_2SITE_3D_FULL, MODEL_NS_CPMG_2SITE_EXPANDED,
MODEL_NS_CPMG_2SITE_STAR, MODEL_NS_CPMG_2SITE_STAR_FULL,
MODEL_NS_R1RHO_2SITE]:
self.inv_relax_time = 1.0 / relax_time
# Special matrices for the multi-quantum CPMG 2-site numerical model.
@@ -823,10 +823,6 @@
k_BA = pA * kex
k_AB = pB * kex
- # This is a vector that contains the initial magnetizations
corresponding to the A and B state transverse magnetizations.
- self.M0[0] = pA
- self.M0[1] = pB
-
# Initialise.
chi2_sum = 0.0
@@ -842,7 +838,7 @@
dwH_frq = dwH[spin_index] * self.frqs[spin_index, frq_index]
# Back calculate the R2eff values.
- r2eff_mq_cr72(r20=R20[r20_index], pA=pA, pB=pB, dw=dw_frq,
dwH=dwH_frq, kex=kex, k_AB=k_AB, k_BA=k_BA, inv_tcpmg=self.inv_relax_time,
tcp=self.tau_cpmg, back_calc=self.back_calc[spin_index, frq_index],
num_points=self.num_disp_points, power=self.power)
+ r2eff_mq_cr72(r20=R20[r20_index], pA=pA, pB=pB, dw=dw_frq,
dwH=dwH_frq, kex=kex, k_AB=k_AB, k_BA=k_BA, cpmg_frqs=self.cpmg_frqs,
tcp=self.tau_cpmg, back_calc=self.back_calc[spin_index, frq_index],
num_points=self.num_disp_points, power=self.power)
# 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):
r21127 - in /branches/relax_disp: lib/dispersion/__init__.py lib/dispersion/mq_cr72.py target_functions/relax_disp.py