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 Meiboom (1961) 2-site on-resonance skewed population R1rho U{M61 skew<http://wiki.nmr-relax.com/M61_skew>} model.
26
27 Description
28 ===========
29
30 This module is for the function, gradient and Hessian of the U{M61 skew<http://wiki.nmr-relax.com/M61_skew>} model.
31
32
33 References
34 ==========
35
36 The model is named after the reference:
37
38 - Meiboom S. (1961). Nuclear magnetic resonance study of the proton transfer in water. I{J. Chem. Phys.}, B{34}, 375-388. (U{DOI: 10.1063/1.1700960<http://dx.doi.org/10.1063/1.1700960>}).
39
40
41 Equations
42 =========
43
44 The equation used is::
45
46 pA^2.pB.delta_omega^2.kex
47 R1rho = R1rho' + -------------------------------------- ,
48 kex^2 + pA^2.delta_omega^2 + omega_1^2
49
50 where R1rho' is the R1rho value in the absence of exchange, kex is the chemical exchange rate constant, pA and pB are the populations of states A and B, delta_omega is the chemical shift difference between the two states, and omega_1 = omega_e is the effective field in the rotating frame.
51
52
53 Links
54 =====
55
56 More information on the M61 skew model can be found in the:
57
58 - U{relax wiki<http://wiki.nmr-relax.com/M61_skew>},
59 - U{relax manual<http://www.nmr-relax.com/manual/The_M61_skew_2_site_fast_exchange_R1_rho_model.html>},
60 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#M61_skew>}.
61 """
62
63
64 from numpy import any, isfinite, min, sum
65 from numpy.ma import fix_invalid, masked_where
66
67
68 -def r1rho_M61b(r1rho_prime=None, pA=None, dw=None, kex=None, spin_lock_fields2=None, back_calc=None):
69 """Calculate the R1rho values for the M61 skew model.
70
71 See the module docstring for details.
72
73
74 @keyword r1rho_prime: The R1rho_prime parameter value (R1rho with no exchange).
75 @type r1rho_prime: numpy float array of rank [NE][NS][NM][NO][ND]
76 @keyword pA: The population of state A.
77 @type pA: float
78 @keyword dw: The chemical exchange difference between states A and B in rad/s.
79 @type dw: numpy float array of rank [NE][NS][NM][NO][ND]
80 @keyword kex: The kex parameter value (the exchange rate in rad/s).
81 @type kex: float
82 @keyword spin_lock_fields2: The R1rho spin-lock field strengths squared (in rad^2.s^-2).
83 @type spin_lock_fields2: numpy float array of rank [NE][NS][NM][NO][ND]
84 @keyword back_calc: The array for holding the back calculated R1rho values. Each element corresponds to the combination of spin lock field.
85 @type back_calc: numpy float array of rank [NE][NS][NM][NO][ND]
86 """
87
88
89 t_numer_zero = False
90 t_denom_zero = False
91
92
93 pB = 1.0 - pA
94
95
96 pA2dw2 = pA**2 * dw**2
97 kex2_pA2dw2 = kex**2 + pA2dw2
98
99
100 numer = pA2dw2 * pB * kex
101
102
103
104 if min(numer) == 0.0:
105 t_numer_zero = True
106 mask_numer_zero = masked_where(numer == 0.0, numer)
107
108
109 denom = kex2_pA2dw2 + spin_lock_fields2
110
111
112
113 mask_denom_zero = denom == 0.0
114 if any(mask_denom_zero):
115 t_denom_zero = True
116 denom[mask_denom_zero] = 1.0
117
118
119 back_calc[:] = r1rho_prime + numer / denom
120
121
122
123 if t_numer_zero:
124 back_calc[mask_numer_zero.mask] = r1rho_prime[mask_numer_zero.mask]
125
126
127 if t_denom_zero:
128 back_calc[mask_denom_zero] = 1e100
129
130
131
132 if not isfinite(sum(back_calc)):
133
134 fix_invalid(back_calc, copy=False, fill_value=1e100)
135