Author: bugman
Date: Fri Oct 11 14:48:41 2013
New Revision: 21068
URL: http://svn.gna.org/viewcvs/relax?rev=21068&view=rev
Log:
Added the 'MQ NS CPMG 2-site' R1rho calculating function to the relax library.
This is the 2-site numeric solution for multi-quantum CPMG-type data.
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.
Added:
branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py
- copied, changed from r21066,
branches/relax_disp/lib/dispersion/ns_cpmg_2site_3d.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=21068&r1=21067&r2=21068&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/__init__.py (original)
+++ branches/relax_disp/lib/dispersion/__init__.py Fri Oct 11 14:48:41 2013
@@ -30,6 +30,7 @@
'lm63_3site',
'm61',
'm61b',
+ 'mq_ns_cpmg_2site',
'ns_cpmg_2site_3d',
'ns_cpmg_2site_expanded',
'ns_cpmg_2site_star',
Copied: branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py (from r21066,
branches/relax_disp/lib/dispersion/ns_cpmg_2site_3d.py)
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py?p2=branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py&p1=branches/relax_disp/lib/dispersion/ns_cpmg_2site_3d.py&r1=21066&r2=21068&rev=21068&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/ns_cpmg_2site_3d.py (original)
+++ branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py Fri Oct 11
14:48:41 2013
@@ -1,6 +1,5 @@
###############################################################################
#
#
-# Copyright (C) 2010-2013 Paul Schanda (https://gna.org/users/pasa)
#
# Copyright (C) 2013 Mathilde Lescanne
#
# Copyright (C) 2013 Dominique Marion
#
# Copyright (C) 2013 Edward d'Auvergne
#
@@ -23,11 +22,11 @@
###############################################################################
# Module docstring.
-"""This function performs a numerical fit of 2-site Bloch-McConnell
equations for CPMG-type experiments.
-
-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. This function was written, initially in MATLAB, in
2010.
-
-This is the model of the numerical solution for the 2-site Bloch-McConnell
equations. It originates as optimization function number 1 from the
fitting_main_kex.py script from Mathilde Lescanne, Paul Schanda, and
Dominique Marion (see
http://thread.gmane.org/gmane.science.nmr.relax.devel/4138,
https://gna.org/task/?7712#comment2 and
https://gna.org/support/download.php?file_id=18262).
+"""This function performs a numerical fit of 2-site Bloch-McConnell
equations for MQ CPMG-type experiments.
+
+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.
+
+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).
"""
# Dependency check module.
@@ -35,7 +34,8 @@
# Python module imports.
from math import fabs, log
-from numpy import dot
+from numpy import array, conj, dot, float64
+from numpy.linalg import matrix_power
if dep_check.scipy_module:
from scipy.linalg import expm
@@ -44,30 +44,86 @@
from lib.float import isNaN
-def r2eff_ns_cpmg_2site_3D(r180x=None, M0=None, r10a=0.0, r10b=0.0,
r20a=None, r20b=None, pA=None, pB=None, dw=None, k_AB=None, k_BA=None,
inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None):
+def populate_m1(matrix=None, r20=None, dw=None, dwH=None, k_AB=None,
k_BA=None):
+ """Matrix for HMQC experiments.
+
+ This corresponds to matrix m1 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 populate_m2(matrix=None, r20=None, dw=None, dwH=None, k_AB=None,
k_BA=None):
+ """Matrix for HMQC experiments.
+
+ This corresponds to matrix m1 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.
This function calculates and stores the R2eff values.
- @keyword r180x: The X-axis pi-pulse propagator.
- @type r180x: numpy float64, rank-2, 7D array
@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 r10a: The R1 value for state A.
- @type r10a: float
- @keyword r10b: The R1 value for state B.
- @type r10b: float
- @keyword r20a: The R2 value for state A in the absence of
exchange.
- @type r20a: float
- @keyword r20b: The R2 value for state B in the absence of
exchange.
- @type r20b: float
+ @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.
@type pA: float
@keyword pB: The population of state B.
@type pB: 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).
@@ -84,25 +140,73 @@
@type power: numpy int16, rank-1 array
"""
- # The matrix R that contains all the contributions to the evolution,
i.e. relaxation, exchange and chemical shift evolution.
- R = rcpmg_3d(R1A=r10a, R1B=r10b, R2A=r20a, R2B=r20b, pA=pA, pB=pB,
dw=dw, k_AB=k_AB, k_BA=k_BA)
+ # Populate the m1 and m2 matrices (only once per function call for
speed).
+ populate_m1(matrix=m1, r20=r20, dw=dw, dwH=dwH, k_AB=k_AB, k_BA=k_BA)
+ populate_m2(matrix=m2, r20=r20, dw=dw, dwH=dwH, k_AB=k_AB, k_BA=k_BA)
# Loop over the time points, back calculating the R2eff values.
for i in range(num_points):
- # Initial magnetisation.
- Mint = M0
-
- # This matrix is a propagator that will evolve the magnetization
with the matrix R for a delay tcp.
- Rexpo = expm(R*tcp[i])
-
- # Loop over the CPMG elements, propagating the magnetisation.
- for j in range(2*power[i]):
- Mint = dot(Rexpo, Mint)
- Mint = dot(r180x, Mint)
- Mint = dot(Rexpo, Mint)
+ # The M1 and M2 matrices.
+ M1 = expm(m1*tcp[i])
+ M2 = expm(m2*tcp[i])
+
+ # The complex conjugates M1* and M2*
+ M1_star = conj(M1)
+ M2_star = conj(M2)
+
+ # 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)
+
+ # Matrices for even n.
+ if power[i] % 2 == 0:
+ # The power factor (only calculate once).
+ fact = power[i] / 2
+
+ # (M1.M2.M2.M1)^(n/2)
+ A = matrix_power(M1_M2_M2_M1, fact)
+
+ # (M2*.M1*.M1*.M2*)^(n/2)
+ B = matrix_power(M2_M1_M1_M2_star, fact)
+
+ # (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.
- Mx = fabs(Mint[1] / pA)
+ 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:
r21068 - in /branches/relax_disp/lib/dispersion: __init__.py mq_ns_cpmg_2site.py