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 math import atan, cos, pi, sin
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 one of the spin-lock fields.
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 for i in range(num_points):
109
110 W = pA*Wa + pB*Wb
111 da = Wa - offset
112 db = Wb - offset
113 d = W - offset
114 waeff2 = spin_lock_fields2[i] + da**2
115 wbeff2 = spin_lock_fields2[i] + db**2
116 weff2 = spin_lock_fields2[i] + d**2
117
118
119 theta = atan(spin_lock_fields[i] / 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 if numer == 0.0:
128 back_calc[i] = R1_cos_theta2 + R1rho_prime_sin_theta2
129 continue
130
131
132 waeff2_wbeff2 = waeff2*wbeff2
133 fact = 1.0 + 2.0*kex2*(pA*waeff2 + pB*wbeff2) / (waeff2_wbeff2 + weff2*kex2)
134 denom = waeff2_wbeff2/weff2 + kex2 - sin_theta2*phi_ex*(fact)
135
136
137 if denom == 0.0:
138 back_calc[i] = 1e100
139 continue
140
141
142 back_calc[i] = R1_cos_theta2 + R1rho_prime_sin_theta2 + sin_theta2 * numer / denom
143