1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
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
59 from math import log
60 from numpy import add, complex, conj, dot
61
62
63 from lib.float import isNaN
64 from lib.linear_algebra.matrix_exponential import matrix_exponential
65 from lib.linear_algebra.matrix_power import square_matrix_power
66
67
68 -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):
69 """The 2-site numerical solution to the Bloch-McConnell equation using complex conjugate matrices.
70
71 This function calculates and stores the R2eff values.
72
73
74 @keyword Rr: The matrix that contains only the R2 relaxation terms ("Redfield relaxation", i.e. non-exchange broadening).
75 @type Rr: numpy complex64, rank-2, 2D array
76 @keyword Rex: The matrix that contains the exchange terms between the two states A and B.
77 @type Rex: numpy complex64, rank-2, 2D array
78 @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).
79 @type RCS: numpy complex64, rank-2, 2D array
80 @keyword R: The matrix that contains all the contributions to the evolution, i.e. relaxation, exchange and chemical shift evolution.
81 @type R: numpy complex64, rank-2, 2D array
82 @keyword M0: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations.
83 @type M0: numpy float64, rank-1, 2D array
84 @keyword r20a: The R2 value for state A in the absence of exchange.
85 @type r20a: float
86 @keyword r20b: The R2 value for state B in the absence of exchange.
87 @type r20b: float
88 @keyword dw: The chemical exchange difference between states A and B in rad/s.
89 @type dw: float
90 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds).
91 @type inv_tcpmg: float
92 @keyword tcp: The tau_CPMG times (1 / 4.nu1).
93 @type tcp: numpy rank-1 float array
94 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies.
95 @type back_calc: numpy rank-1 float array
96 @keyword num_points: The number of points on the dispersion curve, equal to the length of the tcp and back_calc arguments.
97 @type num_points: int
98 @keyword power: The matrix exponential power array.
99 @type power: numpy int16, rank-1 array
100 """
101
102
103 Rr[0, 0] = -r20a
104 Rr[1, 1] = -r20b
105
106
107 RCS[1, 1] = complex(0.0, -dw)
108
109
110 R = add(Rr, Rex)
111 R = add(R, RCS)
112
113
114 cR2 = conj(R) * 2.0
115
116
117 for i in range(num_points):
118
119 eR_tcp = matrix_exponential(R*tcp[i])
120
121
122 prop_2 = dot(dot(eR_tcp, matrix_exponential(cR2*tcp[i])), eR_tcp)
123
124
125 prop_total = square_matrix_power(prop_2, power[i])
126
127
128 Moft = dot(prop_total, M0)
129
130
131 Mx = Moft[0].real / M0[0]
132 if Mx <= 0.0 or isNaN(Mx):
133 back_calc[i] = 1e99
134 else:
135 back_calc[i]= -inv_tcpmg * log(Mx)
136