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 Miloushev and Palmer (2005) 2-site exchange R1rho U{MP05<http://wiki.nmr-relax.com/MP05>} model.
26
27 Description
28 ===========
29
30 This module is for the function, gradient and Hessian of the U{MP05<http://wiki.nmr-relax.com/MP05>} model.
31
32
33 References
34 ==========
35
36 The model is named after the reference:
37
38 - Miloushev V. Z. and Palmer A. (2005). R1rho relaxation for two-site chemical exchange: General approximations and some exact solutions. J. Magn. Reson., 177, 221-227. (U{DOI: 10.1016/j.jmr.2005.07.023<http://dx.doi.org/10.1016/j.jmr.2005.07.023>}).
39
40
41 Code origin
42 ===========
43
44 The code was copied from the U{TP02<http://wiki.nmr-relax.com/TP02>} 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 MP05 model can be found in the:
56
57 - U{relax wiki<http://wiki.nmr-relax.com/MP05>},
58 - U{relax manual<http://www.nmr-relax.com/manual/MP05_2_site_exchange_R1_model.html>},
59 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#MP05>}.
60 """
61
62
63 from numpy import abs, arctan2, array, isfinite, min, sin, sum
64
65
66 -def r1rho_MP05(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 Miloushev and Palmer (2005) 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
103
104 phi_ex = pA * pB * dw**2
105 numer = phi_ex * kex
106
107
108 W = pA*Wa + pB*Wb
109 da = Wa - offset
110 db = Wb - offset
111 d = W - offset
112 waeff2 = spin_lock_fields2 + da**2
113 wbeff2 = spin_lock_fields2 + db**2
114 weff2 = spin_lock_fields2 + d**2
115
116
117 theta = arctan2(spin_lock_fields, d)
118
119
120 sin_theta2 = sin(theta)**2
121 R1_cos_theta2 = R1 * (1.0 - sin_theta2)
122 R1rho_prime_sin_theta2 = r1rho_prime * sin_theta2
123
124
125
126 if numer == 0.0:
127 back_calc[:] = R1_cos_theta2 + R1rho_prime_sin_theta2
128 return
129
130
131 waeff2_wbeff2 = waeff2*wbeff2
132 fact_denom = waeff2_wbeff2 + weff2*kex2
133
134 fact = 1.0 + 2.0*kex2*(pA*waeff2 + pB*wbeff2) / fact_denom
135 denom = waeff2_wbeff2/weff2 + kex2 - sin_theta2*phi_ex*(fact)
136
137
138 R1rho = R1_cos_theta2 + R1rho_prime_sin_theta2 + sin_theta2 * numer / denom
139
140
141
142 if not isfinite(sum(R1rho)):
143 R1rho = array([1e100]*num_points)
144
145 back_calc[:] = R1rho
146