Author: bugman
Date: Wed Jul 17 16:35:25 2013
New Revision: 20361
URL: http://svn.gna.org/viewcvs/relax?rev=20361&view=rev
Log:
Added the 'NS 2-site expanded' R2eff calculating function to the relax
library.
This is the numerical model for the 2-site Bloch-McConnell equations expanded
using Maple by Nikolai
Skrynnikov. It originates as optimization function number 5 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 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_2site_expanded.py
- copied, changed from r20353,
branches/relax_disp/lib/dispersion/ns_2site_star.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=20361&r1=20360&r2=20361&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/__init__.py (original)
+++ branches/relax_disp/lib/dispersion/__init__.py Wed Jul 17 16:35:25 2013
@@ -30,6 +30,7 @@
'm61',
'm61b',
'ns_2site_3d',
+ 'ns_2site_expanded',
'ns_2site_star',
'ns_matrices',
'two_point'
Copied: branches/relax_disp/lib/dispersion/ns_2site_expanded.py (from r20353,
branches/relax_disp/lib/dispersion/ns_2site_star.py)
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/ns_2site_expanded.py?p2=branches/relax_disp/lib/dispersion/ns_2site_expanded.py&p1=branches/relax_disp/lib/dispersion/ns_2site_star.py&r1=20353&r2=20361&rev=20361&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/ns_2site_star.py (original)
+++ branches/relax_disp/lib/dispersion/ns_2site_expanded.py Wed Jul 17
16:35:25 2013
@@ -11,7 +11,7 @@
# (at your option) any later version.
#
#
#
# This program is distributed in the hope that it will be useful,
#
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
#
+# but WITHOUT ANY WARRANTY without even the implied warranty of
#
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#
# GNU General Public License for more details.
#
#
#
@@ -23,9 +23,9 @@
# 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 with
conjugate complex matrices. 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 function is exact, just as the explicit Bloch-McConnell numerical
treatments. It comes from a Maple derivation based on the Bloch-McConnell
equations. It is much faster than the numerical Bloch-McConnell solution.
It was derived by Nikolai Skrynnikov and is provided with his permission.
-The code was submitted at
http://thread.gmane.org/gmane.science.nmr.relax.devel/4132 by Paul Schanda.
+The code originates as optimization function number 5 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).
"""
# Dependency check module.
@@ -33,77 +33,108 @@
# Python module imports.
from math import log
-from numpy import add, complex, conj, dot
+import numpy
+from numpy import add, complex, conj, dot, exp, power, real
from numpy.linalg import matrix_power
if dep_check.scipy_module:
from scipy.linalg import expm
-def r2eff_ns_2site_star(Rr=None, Rex=None, RCS=None, R=None, M0=None,
r20a=None, r20b=None, dw=None, inv_tcpmg=None, tcp=None, back_calc=None,
num_points=None, power=None):
+def r2eff_ns_2site_expanded(r20=None, dw=None, k_AB=None, k_BA=None,
relax_time=None, inv_relax_time=None, tcp=None, back_calc=None,
num_points=None, num_cpmg=None):
"""The 2-site numerical solution to the Bloch-McConnell equation using
complex conjugate matrices.
This function calculates and stores the R2eff values.
- @keyword Rr: The matrix that contains only the R2 relaxation
terms ("Redfield relaxation", i.e. non-exchange broadening).
- @type Rr: numpy complex64, rank-2, 2D array
- @keyword Rex: The matrix that contains the exchange terms
between the two states A and B.
- @type Rex: numpy complex64, rank-2, 2D array
- @keyword RCS: The matrix that contains the chemical shift
evolution. It works here only with X magnetization, and the complex notation
allows to evolve in the transverse plane (x, y).
- @type RCS: numpy complex64, rank-2, 2D array
- @keyword R: The matrix that contains all the contributions
to the evolution, i.e. relaxation, exchange and chemical shift evolution.
- @type R: numpy complex64, rank-2, 2D 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, 2D array
- @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 dw: The chemical exchange difference between states
A and B in rad/s.
- @type dw: 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 r20: The R2 value for both states A and B 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 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 relax_time: The total relaxation time period (in
seconds).
+ @type relax_time: float
+ @keyword inv_relax_time: The inverse of the total relaxation time
period (in inverse seconds).
+ @type inv_relax_time: 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 num_cpmg: The array of numbers of CPMG blocks.
+ @type num_cpmg: numpy int16, rank-1 array
"""
- # The matrix that contains only the R2 relaxation terms ("Redfield
relaxation", i.e. non-exchange broadening).
- Rr[0, 0] = -r20a
- Rr[1, 1] = -r20b
+ # The expansion factors (in numpy array form for all dispersion points).
+ t3 = complex(0, 1)
+ t4 = t3 * dw
+ t5 = k_BA * k_BA
+ t8 = 2.0 * t3 * k_BA * dw
+ t10 = 2.0 * k_BA * k_AB
+ t11 = dw * dw
+ t14 = 2.0 * t3 * k_AB*dw
+ t15 = k_AB * k_AB
+ t17 = cmath.sqrt(t5 - t8 + t10 - t11 + t14 + t15)
+ t21 = exp((-k_BA + t4 - k_AB + t17) * tcp/4.0)
+ t22 = 1.0/t17
+ t28 = exp((-k_BA + t4 - k_AB - t17) * tcp/4.0)
+ t31 = t21 * t22 * k_AB - t28 * t22 * k_AB
+ t33 = cmath.sqrt(t5 + t8 + t10 - t11 - t14 + t15)
+ t34 = k_BA + t4 - k_AB + t33
+ t37 = exp((-k_BA - t4 - k_AB + t33) * tcp/2.0)
+ t39 = 1.0/t33
+ t41 = k_BA + t4 - k_AB - t33
+ t44 = exp((-k_BA - t4 - k_AB - t33) * tcp/2.0)
+ t47 = t34 * t37 * t39/2.0 - t41 * t44 * t39/2.0
+ t49 = k_BA - t4 - k_AB - t17
+ t51 = t21 * t49 * t22
+ t52 = k_BA - t4 - k_AB + t17
+ t54 = t28 * t52 * t22
+ t55 = -t51 + t54
+ t60 = t37 * t39 * k_AB - t44 * t39 * k_AB
+ t62 = t31 * t47 + t55 * t60/2.0
+ t63 = 1.0/k_AB
+ t68 = -t52 * t63 * t51/2.0 + t49 * t63 * t54/2.0
+ t69 = t62 * t68/2.0
+ t72 = t37 * t41 * t39
+ t76 = t44 * t34 * t39
+ t78 = -t34 * t63 * t72/2.0 + t41 * t63 * t76/2.0
+ t80 = -t72 + t76
+ t82 = t31 * t78/2.0 + t55 * t80/4.0
+ t83 = t82 * t55/2.0
+ t88 = t52 * t21 * t22/2.0 - t49 * t28 * t22/2.0
+ t91 = t88 * t47 + t68 * t60/2.0
+ t92 = t91 * t88
+ t95 = t88 * t78/2.0 + t68 * t80/4.0
+ t96 = t95 * t31
+ t97 = t69 + t83
+ t98 = t97 * t97
+ t99 = t92 + t96
+ t102 = t99 * t99
+ t108 = t62 * t88 + t82 * t31
+ t112 = numpy.sqrt(t98 - 2.0 * t99 * t97 + t102 + 4.0 * (t91 * t68/2.0 +
t95 * t55/2.0) * t108)
+ t113 = t69 + t83 - t92 - t96 - t112
+ t115 = num_cpmg / 2.0
+ t116 = power(t69/2.0 + t83/2.0 + t92/2.0 + t96/2.0 + t112/2.0, t115)
+ t118 = 1.0/t112
+ t120 = t69 + t83 - t92 - t96 + t112
+ t122 = power(t69/2.0 + t83/2.0 + t92/2.0 + t96/2.0 - t112/2.0, t115)
+ t127 = 1.0/t108
+ t139 = 1.0/(k_AB + k_BA) * ((-t113 * t116 * t118/2.0 + t120 * t122 *
t118/2.0) * k_BA + (-t113 * t127 * t116 * t120 * t118/2.0 + t120 * t127 *
t122 * t113 * t118/2.0) * k_AB/2.0)
- # The matrix that contains the chemical shift evolution. It works here
only with X magnetization, and the complex notation allows to evolve in the
transverse plane (x, y). The chemical shift for state A is assumed to be
zero.
- RCS[1, 1] = complex(0.0, -dw)
+ # Calculate the initial and final peak intensities.
+ intensity0 = pg
+ intensity = real(t139) * exp(-relax_time * r20)
- # The matrix R that contains all the contributions to the evolution,
i.e. relaxation, exchange and chemical shift evolution.
- add(Rr, Rex, out=R)
- add(R, RCS, out=R)
+ # The magnetisation vector.
+ Mx = intensity / intensity0
- # This is the complex conjugate of the above. It allows the chemical
shift to run in the other direction, i.e. it is used to evolve the shift
after a 180 deg pulse. The factor of 2 is to minimise the number of
multiplications for the prop_2 matrix calculation.
- cR2 = conj(R) * 2.0
-
- # Loop over the time points, back calculating the R2eff values.
+ # Calculate the R2eff using a two-point approximation, i.e. assuming
that the decay is mono-exponential, and store it for each dispersion point.
for i in range(num_points):
- # This matrix is a propagator that will evolve the magnetization
with the matrix R for a delay tcp.
- expm_R_tcp = expm(R*tcp[i])
-
- # This is the propagator for an element of [delay tcp; 180 deg
pulse; 2 times delay tcp; 180 deg pulse; delay tau], i.e. for 2 times
tau-180-tau.
- prop_2 = dot(dot(expm_R_tcp, expm(cR2*tcp[i])), expm_R_tcp)
-
- # Now create the total propagator that will evolve the magnetization
under the CPMG train, i.e. it applies the above tau-180-tau-tau-180-tau so
many times as required for the CPMG frequency under consideration.
- prop_total = matrix_power(prop_2, power[i])
-
- # Now we apply the above propagator to the initial magnetization
vector - resulting in the magnetization that remains after the full CPMG
pulse train. It is called M of t (t is the time after the CPMG train).
- Moft = dot(prop_total, M0)
-
- # The next lines calculate the R2eff using a two-point
approximation, i.e. assuming that the decay is mono-exponential.
- Mx = Moft[0].real / M0[0]
- if Mx == 0.0:
+ if Mx[i] == 0.0:
back_calc[i] = 1e99
else:
- back_calc[i]= -inv_tcpmg * log(Mx)
+ back_calc[i]= -inv_relax_time * log(Mx[i])
r20361 - in /branches/relax_disp/lib/dispersion: __init__.py ns_2site_expanded.py