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 """The Trott and Palmer (2002) 2-site exchange R1rho U{TP02<http://wiki.nmr-relax.com/TP02>} model.
26
27 Description
28 ===========
29
30 This module is for the function, gradient and Hessian of the U{TP02<http://wiki.nmr-relax.com/TP02>} model.
31
32
33 References
34 ==========
35
36 The model is named after the reference:
37
38 - Trott, O. and Palmer, 3rd, A. G. (2002). R1rho relaxation outside of the fast-exchange limit. I{J. Magn. Reson.}, B{154}(1), 157-160. (U{DOI: 10.1006/jmre.2001.2466<http://dx.doi.org/10.1006/jmre.2001.2466>}).
39
40
41 Code origin
42 ===========
43
44 The code originates as the funTrottPalmer.m Matlab file from the sim_all.tar file attached to task #7712 (U{https://web.archive.org/web/https://gna.org/task/?7712}). This is code from Nikolai Skrynnikov and Martin Tollinger.
45
46 Links to the copyright licensing agreements from all authors are:
47
48 - Nikolai Skrynnikov, U{http://article.gmane.org/gmane.science.nmr.relax.devel/4279},
49 - Martin Tollinger, U{http://article.gmane.org/gmane.science.nmr.relax.devel/4276}.
50
51
52 Links
53 =====
54
55 More information on the TP02 model can be found in the:
56
57 - U{relax wiki<http://wiki.nmr-relax.com/TP02>},
58 - U{relax manual<http://www.nmr-relax.com/manual/TP02_2_site_exchange_R1_model.html>},
59 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#TP02>}.
60 """
61
62
63 from numpy import abs, arctan2, array, isfinite, min, sin, sum
64
65
66 -def r1rho_TP02(r1rho_prime=None, omega=None, offset=None, pA=None, pB=None, dw=None, kex=None, R1=0.0, spin_lock_fields=None, spin_lock_fields2=None, back_calc=None, num_points=None):
67 """Calculate the R1rho' values for the TP02 model.
68
69 See the module docstring for details. This is the Trott and Palmer (2002) equation according to Korzhnev (J. Biomol. NMR (2003), 26, 39-48).
70
71
72 @keyword r1rho_prime: The R1rho_prime parameter value (R1rho with no exchange).
73 @type r1rho_prime: float
74 @keyword omega: The chemical shift for the spin in rad/s.
75 @type omega: float
76 @keyword offset: The spin-lock offsets for the data.
77 @type offset: numpy rank-1 float array
78 @keyword pA: The population of state A.
79 @type pA: float
80 @keyword pB: The population of state B.
81 @type pB: float
82 @keyword dw: The chemical exchange difference between states A and B in rad/s.
83 @type dw: float
84 @keyword kex: The kex parameter value (the exchange rate in rad/s).
85 @type kex: float
86 @keyword R1: The R1 relaxation rate.
87 @type R1: float
88 @keyword spin_lock_fields: The R1rho spin-lock field strengths (in rad.s^-1).
89 @type spin_lock_fields: numpy rank-1 float array
90 @keyword spin_lock_fields2: The R1rho spin-lock field strengths squared (in rad^2.s^-2). This is for speed.
91 @type spin_lock_fields2: numpy rank-1 float array
92 @keyword back_calc: The array for holding the back calculated R1rho values. Each element corresponds to the combination of offset and spin lock field.
93 @type back_calc: numpy rank-1 float array
94 @keyword num_points: The number of points on the dispersion curve, equal to the length of the spin_lock_fields and back_calc arguments.
95 @type num_points: int
96 """
97
98
99 Wa = omega
100 Wb = omega + dw
101 kex2 = kex**2
102 W = pA*Wa + pB*Wb
103 da = Wa - offset
104 db = Wb - offset
105 d = W - offset
106 da2 = da**2
107 db2 = db**2
108 d2 = d**2
109
110
111 numer = pA * pB * dw**2 * kex
112
113
114 waeff2 = spin_lock_fields2 + da2
115 wbeff2 = spin_lock_fields2 + db2
116 weff2 = spin_lock_fields2 + d2
117
118
119 theta = arctan2(spin_lock_fields, d)
120
121
122 sin_theta2 = sin(theta)**2
123 R1_cos_theta2 = R1 * (1.0 - sin_theta2)
124 R1rho_prime_sin_theta2 = r1rho_prime * sin_theta2
125
126
127
128 if numer == 0.0:
129 back_calc[:] = R1_cos_theta2 + R1rho_prime_sin_theta2
130 return
131
132
133 denom = waeff2 * wbeff2 / weff2 + kex2
134
135
136
137 R1rho = R1_cos_theta2 + R1rho_prime_sin_theta2 + sin_theta2 * numer / denom
138
139
140
141 if not isfinite(sum(R1rho)):
142 R1rho = array([1e100]*num_points)
143
144 back_calc[:] = R1rho
145