| Trees | Indices | Help |
|
|---|
|
|
1 ###############################################################################
2 # #
3 # Copyright (C) 2000-2001 Nikolai Skrynnikov #
4 # Copyright (C) 2000-2001 Martin Tollinger #
5 # Copyright (C) 2010-2013 Paul Schanda (https://web.archive.org/web/https://gna.org/users/pasa) #
6 # Copyright (C) 2013 Mathilde Lescanne #
7 # Copyright (C) 2013 Dominique Marion #
8 # Copyright (C) 2013-2014 Edward d'Auvergne #
9 # #
10 # This file is part of the program relax (http://www.nmr-relax.com). #
11 # #
12 # This program is free software: you can redistribute it and/or modify #
13 # it under the terms of the GNU General Public License as published by #
14 # the Free Software Foundation, either version 3 of the License, or #
15 # (at your option) any later version. #
16 # #
17 # This program is distributed in the hope that it will be useful, #
18 # but WITHOUT ANY WARRANTY; without even the implied warranty of #
19 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
20 # GNU General Public License for more details. #
21 # #
22 # You should have received a copy of the GNU General Public License #
23 # along with this program. If not, see <http://www.gnu.org/licenses/>. #
24 # #
25 ###############################################################################
26
27 # Module docstring.
28 """The numerical fit of 2-site Bloch-McConnell equations for CPMG-type experiments, the U{NS CPMG 2-site star<http://wiki.nmr-relax.com/NS_CPMG_2-site_star>} and U{NS CPMG 2-site star full<http://wiki.nmr-relax.com/NS_CPMG_2-site_star_full>} models.
29
30 Description
31 ===========
32
33 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. 2004 I{Chem. Rev.}, B{104}, 3623-3640. This function was written, initially in MATLAB, in 2010.
34
35
36 Code origin
37 ===========
38
39 The code was submitted at U{http://thread.gmane.org/gmane.science.nmr.relax.devel/4132} by Paul Schanda.
40
41
42 Links
43 =====
44
45 More information on the NS CPMG 2-site star model can be found in the:
46
47 - U{relax wiki<http://wiki.nmr-relax.com/NS_CPMG_2-site_star>},
48 - U{relax manual<http://www.nmr-relax.com/manual/reduced_NS_2_site_star_CPMG_model.html>},
49 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_CPMG_2-site_star>}.
50
51 More information on the NS CPMG 2-site star full model can be found in the:
52
53 - U{relax wiki<http://wiki.nmr-relax.com/NS_CPMG_2-site_star_full>},
54 - U{relax manual<http://www.nmr-relax.com/manual/full_NS_2_site_star_CPMG_model.html>},
55 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_CPMG_2-site_star_full>}.
56 """
57
58 # Dependency check module.
59 import dep_check
60
61 # Python module imports.
62 from math import log
63 from numpy import add, complex, conj, dot
64
65 # relax module imports.
66 from lib.float import isNaN
67 from lib.linear_algebra.matrix_exponential import matrix_exponential
68 from lib.linear_algebra.matrix_power import square_matrix_power
69
70
71 -def r2eff_ns_cpmg_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):
72 """The 2-site numerical solution to the Bloch-McConnell equation using complex conjugate matrices.
73
74 This function calculates and stores the R2eff values.
75
76
77 @keyword Rr: The matrix that contains only the R2 relaxation terms ("Redfield relaxation", i.e. non-exchange broadening).
78 @type Rr: numpy complex64, rank-2, 2D array
79 @keyword Rex: The matrix that contains the exchange terms between the two states A and B.
80 @type Rex: numpy complex64, rank-2, 2D array
81 @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).
82 @type RCS: numpy complex64, rank-2, 2D array
83 @keyword R: The matrix that contains all the contributions to the evolution, i.e. relaxation, exchange and chemical shift evolution.
84 @type R: numpy complex64, rank-2, 2D array
85 @keyword M0: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations.
86 @type M0: numpy float64, rank-1, 2D array
87 @keyword r20a: The R2 value for state A in the absence of exchange.
88 @type r20a: float
89 @keyword r20b: The R2 value for state B in the absence of exchange.
90 @type r20b: float
91 @keyword dw: The chemical exchange difference between states A and B in rad/s.
92 @type dw: float
93 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds).
94 @type inv_tcpmg: float
95 @keyword tcp: The tau_CPMG times (1 / 4.nu1).
96 @type tcp: numpy rank-1 float array
97 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies.
98 @type back_calc: numpy rank-1 float array
99 @keyword num_points: The number of points on the dispersion curve, equal to the length of the tcp and back_calc arguments.
100 @type num_points: int
101 @keyword power: The matrix exponential power array.
102 @type power: numpy int16, rank-1 array
103 """
104
105 # The matrix that contains only the R2 relaxation terms ("Redfield relaxation", i.e. non-exchange broadening).
106 Rr[0, 0] = -r20a
107 Rr[1, 1] = -r20b
108
109 # 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.
110 RCS[1, 1] = complex(0.0, -dw)
111
112 # The matrix R that contains all the contributions to the evolution, i.e. relaxation, exchange and chemical shift evolution.
113 R = add(Rr, Rex)
114 R = add(R, RCS)
115
116 # 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.
117 cR2 = conj(R) * 2.0
118
119 # Loop over the time points, back calculating the R2eff values.
120 for i in range(num_points):
121 # This matrix is a propagator that will evolve the magnetization with the matrix R for a delay tcp.
122 eR_tcp = matrix_exponential(R*tcp[i])
123
124 # 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.
125 prop_2 = dot(dot(eR_tcp, matrix_exponential(cR2*tcp[i])), eR_tcp)
126
127 # 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.
128 prop_total = square_matrix_power(prop_2, power[i])
129
130 # 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).
131 Moft = dot(prop_total, M0)
132
133 # The next lines calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential.
134 Mx = Moft[0].real / M0[0]
135 if Mx <= 0.0 or isNaN(Mx):
136 back_calc[i] = 1e99
137 else:
138 back_calc[i]= -inv_tcpmg * log(Mx)
139
| Trees | Indices | Help |
|
|---|
| Generated by Epydoc 3.0.1 on Mon Mar 17 15:12:25 2014 | http://epydoc.sourceforge.net |