Author: bugman
Date: Mon Dec 9 17:30:16 2013
New Revision: 21916
URL: http://svn.gna.org/viewcvs/relax?rev=21916&view=rev
Log:
Modified the 'NS R1rho 2-site' dispersion model to match the 'NS R1rho
3-site' models.
The 6D evolution matrix indices have been rearranged to match the 9D matrix
indices. The tilt angle
for the initial magnetisation is no longer that for the average offset but
that of state A, as was
changed for the 'NS R1rho 3-site' models earlier (r21911). The system test
was therefore updated
for the slightly different behaviour.
Modified:
trunk/docs/latex/dispersion.tex
trunk/lib/dispersion/ns_matrices.py
trunk/lib/dispersion/ns_r1rho_2site.py
trunk/target_functions/relax_disp.py
trunk/test_suite/system_tests/relax_disp.py
Modified: trunk/docs/latex/dispersion.tex
URL:
http://svn.gna.org/viewcvs/relax/trunk/docs/latex/dispersion.tex?rev=21916&r1=21915&r2=21916&view=diff
==============================================================================
--- trunk/docs/latex/dispersion.tex (original)
+++ trunk/docs/latex/dispersion.tex Mon Dec 9 17:30:16 2013
@@ -1111,20 +1111,20 @@
where
\begin{align}
- M_0 &= \begin{pmatrix} \sin{\theta} \\ 0 \\ 0 \\ 0 \\ \cos{\theta} \\
0 \end{pmatrix}, \\
- \theta &= \arctan \left( \frac{\omegaone}{\aveoffset} \right).
+ M_0 &= \begin{pmatrix} \sin{\theta} \\ 0 \\ \cos{\theta} \\ 0 \\ 0 \\
0 \end{pmatrix}, \\
+ \theta &= \arctan \left( \frac{\omegaone}{\offsetA} \right).
\end{align}
The relaxation evolution matrix is defined as
\begin{equation}
- R = -\begin{pmatrix}
- \Ronerhoprime+\kAB & -\kBA & \delta_A & 0
& 0 & 0 \\
- -\kAB & \Ronerhoprime+\kBA & 0 &
\delta_B & 0 & 0 \\
- -\delta_A & 0 & \Ronerhoprime+\kAB &
-\kBA & \omegaone & 0 \\
- 0 & -\delta_B & -\kAB &
\Ronerhoprime+\kBA & 0 & \omegaone \\
- 0 & 0 & -\omegaone & 0
& \Rone+\kAB & -\kBA \\
- 0 & 0 & 0 &
-\omegaone & -\kAB & \Rone+\kBA \\
- \end{pmatrix},
+ R = \begin{pmatrix}
+ -\Ronerhoprime-\kAB & -\delta_A & 0 & \kBA
& 0 & 0 \\
+ \delta_A & -\Ronerhoprime-\kAB & -\omegaone & 0
& \kBA & 0 \\
+ 0 & \omegaone & -\Rone-\kAB & 0
& 0 & \kBA \\
+ \kAB & 0 & 0 &
-\Ronerhoprime-\kBA & -\delta_B & 0 \\
+ 0 & \kAB & 0 & \delta_B
& -\Ronerhoprime-\kBA & -\omegaone \\
+ 0 & 0 & \kAB & 0
& \omegaone & -\Rone-\kBA \\
+ \end{pmatrix},
\end{equation}
where $\delta_{A,B}$ is defined in Equations~\ref{eq: deltaA} and~\ref{eq:
deltaB}.
Modified: trunk/lib/dispersion/ns_matrices.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/lib/dispersion/ns_matrices.py?rev=21916&r1=21915&r2=21916&view=diff
==============================================================================
--- trunk/lib/dispersion/ns_matrices.py (original)
+++ trunk/lib/dispersion/ns_matrices.py Mon Dec 9 17:30:16 2013
@@ -249,43 +249,56 @@
matrix[8, 5] = k_BC
-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):
+def rr1rho_3d(matrix=None, R1=None, r1rho_prime=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
+ @keyword matrix: The matrix to fill.
+ @type matrix: numpy rank-2 6D array
+ @keyword R1: The longitudinal, spin-lattice relaxation rate.
+ @type R1: float
+ @keyword r1rho_prime: The R1rho transverse, spin-spin relaxation rate
in the absence of exchange.
+ @type r1rho_prime: 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
+ """
+
+ # The AB auto-block.
+ matrix[0, 0] = -r1rho_prime - k_AB
+ matrix[0, 1] = -wA
+ matrix[1, 0] = wA
+ matrix[1, 1] = -r1rho_prime - k_AB
+ matrix[1, 2] = -w1
+ matrix[2, 1] = w1
+ matrix[2, 2] = -R1 - k_AB
+
+ # The BA auto-block.
+ matrix[3, 3] = -r1rho_prime - k_BA
+ matrix[3, 4] = -wB
+ matrix[4, 3] = wB
+ matrix[4, 4] = -r1rho_prime - k_BA
+ matrix[4, 5] = -w1
+ matrix[5, 4] = w1
+ matrix[5, 5] = -R1 - k_BA
+
+ # The AB cross-block.
+ matrix[0, 3] = k_BA
+ matrix[1, 4] = k_BA
+ matrix[2, 5] = k_BA
+ matrix[3, 0] = k_AB
+ matrix[4, 1] = k_AB
+ matrix[5, 2] = k_AB
Modified: trunk/lib/dispersion/ns_r1rho_2site.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/lib/dispersion/ns_r1rho_2site.py?rev=21916&r1=21915&r2=21916&view=diff
==============================================================================
--- trunk/lib/dispersion/ns_r1rho_2site.py (original)
+++ trunk/lib/dispersion/ns_r1rho_2site.py Mon Dec 9 17:30:16 2013
@@ -46,7 +46,7 @@
from lib.linear_algebra.matrix_exponential import matrix_exponential
-def ns_r1rho_2site(M0=None, r1rho_prime=None, omega=None, offset=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):
+def ns_r1rho_2site(M0=None, matrix=None, r1rho_prime=None, omega=None,
offset=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 R1rho values.
@@ -54,6 +54,8 @@
@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 matrix: A numpy array to be populated to create the
evolution matrix.
+ @type matrix: numpy rank-2, 6D float64 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.
@@ -94,16 +96,16 @@
# Loop over the time points, back calculating the R2eff values.
for i in range(num_points):
- # 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=spin_lock_fields[i], k_AB=k_AB, k_BA=k_BA)
+ # The matrix that contains all the contributions to the evolution,
i.e. relaxation, exchange and chemical shift evolution.
+ rr1rho_3d(matrix=matrix, R1=r1, r1rho_prime=r1rho_prime, pA=pA,
pB=pB, wA=dA, wB=dB, w1=spin_lock_fields[i], k_AB=k_AB, k_BA=k_BA)
# The following lines rotate the magnetization previous to spin-lock
into the weff frame.
- theta = atan(spin_lock_fields[i]/d)
- M0[0] = sin(theta)
- M0[4] = cos(theta)
+ theta = atan(spin_lock_fields[i]/dA)
+ M0[0] = sin(theta) # The A state initial X magnetisation.
+ M0[2] = cos(theta) # The A state initial Z magnetisation.
# This matrix is a propagator that will evolve the magnetization
with the matrix R.
- Rexpo = matrix_exponential(R*relax_time)
+ Rexpo = matrix_exponential(matrix*relax_time)
# Magnetization evolution.
MA = dot(M0, dot(Rexpo, M0))
Modified: trunk/target_functions/relax_disp.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/target_functions/relax_disp.py?rev=21916&r1=21915&r2=21916&view=diff
==============================================================================
--- trunk/target_functions/relax_disp.py (original)
+++ trunk/target_functions/relax_disp.py Mon Dec 9 17:30:16 2013
@@ -328,6 +328,8 @@
elif model in [MODEL_NS_MMQ_3SITE, MODEL_NS_MMQ_3SITE_LINEAR]:
self.m1 = zeros((3, 3), complex64)
self.m2 = zeros((3, 3), complex64)
+ elif model == MODEL_NS_R1RHO_2SITE:
+ self.matrix = zeros((6, 6), float64)
elif model in [MODEL_NS_R1RHO_3SITE, MODEL_NS_R1RHO_3SITE_LINEAR]:
self.matrix = zeros((9, 9), float64)
@@ -1596,7 +1598,7 @@
# Loop over the offsets.
for oi in range(self.num_offsets[0][si][mi]):
# Back calculate the R2eff values.
- ns_r1rho_2site(M0=self.M0,
r1rho_prime=r1rho_prime[r20_index], omega=self.chemical_shifts[0][si][mi],
offset=self.offset[0][si][mi][oi], r1=self.r1[si, mi], pA=pA, pB=pB,
dw=dw_frq, k_AB=k_AB, k_BA=k_BA,
spin_lock_fields=self.spin_lock_omega1[0][mi][oi],
relax_time=self.relax_times[0][mi],
inv_relax_time=self.inv_relax_times[0][mi],
back_calc=self.back_calc[0][si][mi][oi],
num_points=self.num_disp_points[0][si][mi][oi])
+ ns_r1rho_2site(M0=self.M0, matrix=self.matrix,
r1rho_prime=r1rho_prime[r20_index], omega=self.chemical_shifts[0][si][mi],
offset=self.offset[0][si][mi][oi], r1=self.r1[si, mi], pA=pA, pB=pB,
dw=dw_frq, k_AB=k_AB, k_BA=k_BA,
spin_lock_fields=self.spin_lock_omega1[0][mi][oi],
relax_time=self.relax_times[0][mi],
inv_relax_time=self.inv_relax_times[0][mi],
back_calc=self.back_calc[0][si][mi][oi],
num_points=self.num_disp_points[0][si][mi][oi])
# 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 di in range(self.num_disp_points[0][si][mi][oi]):
Modified: trunk/test_suite/system_tests/relax_disp.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/test_suite/system_tests/relax_disp.py?rev=21916&r1=21915&r2=21916&view=diff
==============================================================================
--- trunk/test_suite/system_tests/relax_disp.py (original)
+++ trunk/test_suite/system_tests/relax_disp.py Mon Dec 9 17:30:16 2013
@@ -2987,20 +2987,20 @@
print("%-20s %20.15g %20.15g\n" % ("chi2", spin1.chi2, spin2.chi2))
# Checks for residue :1.
- self.assertAlmostEqual(spin1.r2[r20_key1], 8.50487445281754, 4)
- self.assertAlmostEqual(spin1.r2[r20_key2], 13.4707451492195, 4)
- self.assertAlmostEqual(spin1.pA, 0.864789608883791, 4)
- self.assertAlmostEqual(spin1.dw, 8.86964329737465, 4)
- self.assertAlmostEqual(spin1.kex/1000, 1195.60566701785/1000, 4)
- self.assertAlmostEqual(spin1.chi2, 3.00587274614435, 4)
+ self.assertAlmostEqual(spin1.r2[r20_key1], 8.50207717367548, 4)
+ self.assertAlmostEqual(spin1.r2[r20_key2], 13.4680429589972, 4)
+ self.assertAlmostEqual(spin1.pA, 0.864523128842393, 4)
+ self.assertAlmostEqual(spin1.dw, 8.85204828994151, 4)
+ self.assertAlmostEqual(spin1.kex/1000, 1199.56359549637/1000, 4)
+ self.assertAlmostEqual(spin1.chi2, 2.99182130153514, 4)
# Checks for residue :2.
- self.assertAlmostEqual(spin2.r2[r20_key1], 10.2880464246005, 4)
- self.assertAlmostEqual(spin2.r2[r20_key2], 16.2875959092412, 4)
- self.assertAlmostEqual(spin2.pA, 0.843381553940209, 4)
- self.assertAlmostEqual(spin2.dw, 9.59767243628999, 4)
- self.assertAlmostEqual(spin2.kex/1000, 1491.81723238519/1000, 4)
- self.assertAlmostEqual(spin2.chi2, 0.000861055839088852, 4)
+ self.assertAlmostEqual(spin2.r2[r20_key1], 10.1334196530849, 4)
+ self.assertAlmostEqual(spin2.r2[r20_key2], 16.140167863407, 4)
+ self.assertAlmostEqual(spin2.pA, 0.829988381197468, 4)
+ self.assertAlmostEqual(spin2.dw, 9.5657894936005, 4)
+ self.assertAlmostEqual(spin2.kex/1000, 1404.76852145933/1000, 4)
+ self.assertAlmostEqual(spin2.chi2, 0.000150052743893402, 4)
def test_tp02_data_to_mp05(self):
r21916 - in /trunk: docs/latex/ lib/dispersion/ target_functions/ test_suite/system_tests/