Author: bugman
Date: Fri Aug 30 11:29:38 2013
New Revision: 20723
URL: http://svn.gna.org/viewcvs/relax?rev=20723&view=rev
Log:
Added the 'NS R1rho 2-site' R1rho calculating function to the relax library.
This is the numerical model for the 2-site Bloch-McConnell equations for
R1rho data.
This code originates from the Skrynikov & Tollinger code (the sim_all.tar file
https://gna.org/support/download.php?file_id=18404 attached to
https://gna.org/task/?7712#comment5).
Specifically the funNumrho.m file.
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/ns_r1rho_2site.py
- copied, changed from r20720,
branches/relax_disp/lib/dispersion/ns_cpmg_2site_3d.py
Modified:
branches/relax_disp/lib/dispersion/__init__.py
branches/relax_disp/lib/dispersion/ns_matrices.py
Modified: branches/relax_disp/lib/dispersion/__init__.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/__init__.py?rev=20723&r1=20722&r2=20723&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/__init__.py (original)
+++ branches/relax_disp/lib/dispersion/__init__.py Fri Aug 30 11:29:38 2013
@@ -34,6 +34,7 @@
'ns_cpmg_2site_3d',
'ns_cpmg_2site_expanded',
'ns_cpmg_2site_star',
+ 'ns_r1rho_2site',
'ns_matrices',
'tp02',
'two_point'
Modified: branches/relax_disp/lib/dispersion/ns_matrices.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/ns_matrices.py?rev=20723&r1=20722&r2=20723&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/ns_matrices.py (original)
+++ branches/relax_disp/lib/dispersion/ns_matrices.py Fri Aug 30 11:29:38 2013
@@ -1,5 +1,7 @@
###############################################################################
#
#
+# Copyright (C) 2000-2001 Nikolai Skrynnikov
#
+# Copyright (C) 2000-2001 Martin Tollinger
#
# Copyright (C) 2013 Mathilde Lescanne
#
# Copyright (C) 2013 Dominique Marion
#
# Copyright (C) 2013 Edward d'Auvergne
#
@@ -156,3 +158,45 @@
# Return the matrix.
return temp
+
+
+def rr1rho_3d(R1A=None, R1=None, Rinf=None, pA=None, pB=None, wA=None,
wB=None, w1=None, k_AB=None, k_BA=None):
+ """Definition of the 3D exchange matrix.
+
+ This code originates from the funNumrho.m file from the Skrynikov &
Tollinger code (the sim_all.tar file
https://gna.org/support/download.php?file_id=18404 attached to
https://gna.org/task/?7712#comment5).
+
+
+ @keyword R1: The longitudinal, spin-lattice relaxation rate.
+ @type R1: float
+ @keyword Rinf: The R1rho transverse, spin-spin relaxation rate in the
absence of exchange.
+ @type Rinf: float
+ @keyword pA: The population of state A.
+ @type pA: float
+ @keyword pB: The population of state B.
+ @type pB: float
+ @keyword wA: The chemical shift offset of state A from the spin-lock.
+ @type wA: float
+ @keyword wB: The chemical shift offset of state A from the spin-lock.
+ @type wB: float
+ @keyword w1: The spin-lock field strength in rad/s.
+ @type w1: float
+ @keyword k_AB: The forward exchange rate from state A to state B.
+ @type k_AB: float
+ @keyword k_BA: The reverse exchange rate from state B to state A.
+ @type k_BA: float
+ @return: The relaxation matrix.
+ @rtype: numpy rank-2, 7D array
+ """
+
+ # Create the matrix.
+ temp = -matrix([
+ [ Rinf+k_AB, -k_BA, wA, 0.0, 0.0, 0.0],
+ [ -k_AB, Rinf+k_BA, 0.0, wB, 0.0, 0.0],
+ [ -wA, 0.0, Rinf+k_AB, -k_BA, w1, 0.0],
+ [ 0.0, -wB, -k_AB, Rinf+k_BA, 0.0, w1],
+ [ 0.0, 0.0, -w1, 0.0, R1+k_AB, -k_BA],
+ [ 0.0, 0.0, 0.0, -w1, -k_AB, R1+k_BA]
+ ])
+
+ # Return the matrix.
+ return temp
Copied: branches/relax_disp/lib/dispersion/ns_r1rho_2site.py (from r20720,
branches/relax_disp/lib/dispersion/ns_cpmg_2site_3d.py)
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/ns_r1rho_2site.py?p2=branches/relax_disp/lib/dispersion/ns_r1rho_2site.py&p1=branches/relax_disp/lib/dispersion/ns_cpmg_2site_3d.py&r1=20720&r2=20723&rev=20723&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/ns_cpmg_2site_3d.py (original)
+++ branches/relax_disp/lib/dispersion/ns_r1rho_2site.py Fri Aug 30 11:29:38
2013
@@ -1,5 +1,7 @@
###############################################################################
#
#
+# Copyright (C) 2000-2001 Nikolai Skrynnikov
#
+# Copyright (C) 2000-2001 Martin Tollinger
#
# Copyright (C) 2010-2013 Paul Schanda (https://gna.org/users/pasa)
#
# Copyright (C) 2013 Mathilde Lescanne
#
# Copyright (C) 2013 Dominique Marion
#
@@ -23,87 +25,99 @@
###############################################################################
# Module docstring.
-"""This function performs a numerical fit of 2-site Bloch-McConnell
equations for CPMG-type experiments.
+"""This function performs a numerical fit of 2-site Bloch-McConnell
equations for R1rho-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 is the model of the numerical solution for the 2-site Bloch-McConnell
equations. It originates from the funNumrho.m file from the Skrynikov &
Tollinger code (the sim_all.tar file
https://gna.org/support/download.php?file_id=18404 attached to
https://gna.org/task/?7712#comment5).
"""
# Dependency check module.
import dep_check
# Python module imports.
-from math import fabs, log
+from math import atan, fabs, log
from numpy import dot
if dep_check.scipy_module:
from scipy.linalg import expm
# relax module imports.
-from lib.dispersion.ns_matrices import rcpmg_3d
+from lib.dispersion.ns_matrices import rr1rho_3d
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):
- """The 2-site numerical solution to the Bloch-McConnell equation.
+def r2eff_ns_cpmg_2site_3D(M0=None, r1rho_prime=None, omega=None, r1=0.0,
pA=None, pB=None, dw=None, k_AB=None, k_BA=None, spin_lock_fields=None,
relax_time=None, inv_relax_time=None, back_calc=None, num_points=None):
+ """The 2-site numerical solution to the Bloch-McConnell equation for
R1rho data.
- This function calculates and stores the R2eff values.
+ This function calculates and stores the R1rho 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 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 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 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.
- @type back_calc: numpy rank-1 float array
- @keyword num_points: The number of points on the dispersion curve,
equal to the length of the tcp and back_calc arguments.
- @type num_points: int
- @keyword power: The matrix exponential power array.
- @type power: numpy int16, rank-1 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 r1rho_prime: The R1rho_prime parameter value (R1rho with
no exchange).
+ @type r1rho_prime: float
+ @keyword omega: The chemical shift for the spin in rad/s.
+ @type omega: float
+ @keyword r1: The R1 relaxation rate.
+ @type r1: 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 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 spin_lock_fields: The R1rho spin-lock field strengths in Hz.
+ @type spin_lock_fields: numpy rank-1 float array
+ @keyword relax_time: The total relaxation time period for each
spin-lock field strength (in seconds).
+ @type relax_time: numpy rank-1 float array
+ @keyword inv_relax_time: The inverse of the relaxation time period
for each spin-lock field strength (in inverse seconds). This is used for
faster calculations.
+ @type inv_relax_time: 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.
+ @type back_calc: numpy rank-1 float array
+ @keyword num_points: The number of points on the dispersion
curve, equal to the length of the tcp and back_calc arguments.
+ @type num_points: int
"""
- # 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)
+ # Repetitive calculations (to speed up calculations).
+ Wa = omega # Larmor frequency [s^-1].
+ Wb = omega + dw # Larmor frequency [s^-1].
+ kex2 = kex**2
# Loop over the time points, back calculating the R2eff values.
for i in range(num_points):
- # Initial magnetisation.
- Mint = M0
+ Wsl = offset[i] # Larmor frequency of spin lock
[s^-1].
+ w1 = spin_lock_fields[i] * 2.0 * pi # Spin-lock field strength.
+ dA = Wa - Wsl # Offset of spin-lock from A.
+ dB = Wb - Wsl # Offset of spin-lock from B.
- # This matrix is a propagator that will evolve the magnetization
with the matrix R for a delay tcp.
- Rexpo = expm(R*tcp[i])
+ # The matrix R that contains all the contributions to the evolution,
i.e. relaxation, exchange and chemical shift evolution.
+ R = rr1rho_3d(R1=r1, Rinf=r1rho_prime, pA=pA, pB=pB, wA=dA, wB=dB,
w1=w1, k_AB=k_AB, k_BA=k_BA)
- # 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 following lines rotate the magnetization previous to spin-lock
into the weff frame.
+ # A new M0 is obtained: M0 = [sin(thetaA)*pA sin(thetaB)*pB 0 0
cos(thetaA)*pA cos(thetaB)*pB].
+ thetaA = atan(w1/dA)
+ thetaB = atan(w1/dB)
+ M0[0] = sin(thetaA) * pA
+ M0[1] = sin(thetaB) * pB
+ M0[4] = cos(thetaA) * pA
+ M0[5] = cos(thetaB) * pB
- # 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)
- if Mx <= 0.0 or isNaN(Mx):
+ # This matrix is a propagator that will evolve the magnetization
with the matrix R.
+ Rexpo = expm(R*relax_time[i])
+
+ # Magnetization evolution.
+ Moft = dot(Rexpo, M0)
+ MAx = Moft[0].real / pA
+ MAy = Moft[2].real / pA
+ MAz = Moft[4].real / pA
+ MA = sqrt(MAx**2 + MAy**2 + MAz**2) # For spin A, is equal to 1
if nothing happens (no relaxation).
+
+ # The next lines calculate the R1rho using a two-point
approximation, i.e. assuming that the decay is mono-exponential.
+ if MA <= 0.0 or isNaN(MA):
back_calc[i] = 1e99
else:
- back_calc[i]= -inv_tcpmg * log(Mx)
+ back_calc[i]= -inv_relax_time[i] * log(MA)
+
+
r20723 - in /branches/relax_disp/lib/dispersion: __init__.py ns_matrices.py ns_r1rho_2site.py