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, Abergel and Palmer (2003) off-resonance 2-site exchange R1rho U{TAP03<http://wiki.nmr-relax.com/TAP03>} model.
26
27 Description
28 ===========
29
30 This module is for the function, gradient and Hessian of the U{TAP03<http://wiki.nmr-relax.com/TAP03>} model.
31
32
33 References
34 ==========
35
36 The model is named after the reference:
37
38 - Trott, O., Abergel, D., and Palmer, A. (2003). An average-magnetization analysis of R1rho relaxation outside of the fast exchange. I{Mol. Phys.}, B{101}, 753-763. (U{DOI: 10.1080/0026897021000054826<http://dx.doi.org/10.1080/0026897021000054826>}).
39
40
41 Code origin
42 ===========
43
44 The code was copied from the U{TP02<http://wiki.nmr-relax.com/TP02>} model (via the U{MP05<http://wiki.nmr-relax.com/MP05>} model), hence it 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 TAP03 model can be found in the:
56
57 - U{relax wiki<http://wiki.nmr-relax.com/TAP03>},
58 - U{relax manual<http://www.nmr-relax.com/manual/TAP03_2_site_exchange_R1_model.html>},
59 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#TAP03>}.
60 """
61
62
63 from numpy import arctan2, array, isfinite, min, sin, sqrt, sum
64
65
66 -def r1rho_TAP03(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, Abergel and Palmer (2003) equation.
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
104
105 phi_ex = pA * pB * dw**2
106 numer = phi_ex * kex
107
108
109 da = Wa - offset
110 db = Wb - offset
111 d = W - offset
112
113
114 sigma = pB*da + pA*db
115 sigma2 = sigma**2
116
117 gamma = 1.0 + phi_ex*(sigma2 - kex2 + spin_lock_fields2) / (sigma2 + kex2 + spin_lock_fields2)**2
118
119
120 if min(gamma) < 0.0:
121 back_calc[:] = array([1e100]*num_points)
122 return
123
124
125 waeff2 = gamma*spin_lock_fields2 + da**2
126 wbeff2 = gamma*spin_lock_fields2 + db**2
127 weff2 = gamma*spin_lock_fields2 + d**2
128
129
130 theta = arctan2(spin_lock_fields, d)
131 hat_theta = arctan2(sqrt(gamma)*spin_lock_fields, d)
132
133
134 sin_theta2 = sin(theta)**2
135 hat_sin_theta2 = sin(hat_theta)**2
136 R1_cos_theta2 = R1 * (1.0 - sin_theta2)
137 R1rho_prime_sin_theta2 = r1rho_prime * sin_theta2
138
139
140
141 if numer == 0.0:
142 back_calc[:] = R1_cos_theta2 + R1rho_prime_sin_theta2
143 return
144
145
146 denom = waeff2*wbeff2/weff2 + kex2 - 2.0*hat_sin_theta2*phi_ex + (1.0 - gamma)*spin_lock_fields2
147
148
149 R1rho = R1_cos_theta2 + R1rho_prime_sin_theta2 + hat_sin_theta2 * numer / denom / gamma
150
151
152
153 if not isfinite(sum(R1rho)):
154 R1rho = array([1e100]*num_points)
155
156 back_calc[:] = R1rho
157