Package lib :: Package dispersion :: Module ns_cpmg_2site_3d
[hide private]
[frames] | no frames]

Source Code for Module lib.dispersion.ns_cpmg_2site_3d

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2010-2013 Paul Schanda                                        # 
  4  # Copyright (C) 2013 Mathilde Lescanne                                        # 
  5  # Copyright (C) 2013 Dominique Marion                                         # 
  6  # Copyright (C) 2013-2014 Edward d'Auvergne                                   # 
  7  # Copyright (C) 2014 Troels E. Linnet                                         # 
  8  #                                                                             # 
  9  # This file is part of the program relax (http://www.nmr-relax.com).          # 
 10  #                                                                             # 
 11  # This program is free software: you can redistribute it and/or modify        # 
 12  # it under the terms of the GNU General Public License as published by        # 
 13  # the Free Software Foundation, either version 3 of the License, or           # 
 14  # (at your option) any later version.                                         # 
 15  #                                                                             # 
 16  # This program is distributed in the hope that it will be useful,             # 
 17  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 18  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 19  # GNU General Public License for more details.                                # 
 20  #                                                                             # 
 21  # You should have received a copy of the GNU General Public License           # 
 22  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 23  #                                                                             # 
 24  ############################################################################### 
 25   
 26  # Module docstring. 
 27  """The numerical fit of 2-site Bloch-McConnell equations for CPMG-type experiments, the U{NS CPMG 2-site 3D<http://wiki.nmr-relax.com/NS_CPMG_2-site_3D>} and U{NS CPMG 2-site 3D full<http://wiki.nmr-relax.com/NS_CPMG_2-site_3D_full>} models. 
 28   
 29  Description 
 30  =========== 
 31   
 32  The function uses an explicit matrix that contains relaxation, exchange and chemical shift terms.  It does the 180deg pulses in the CPMG train.  The approach of Bloch-McConnell can be found in chapter 3.1 of Palmer, A. G. Chem Rev 2004, 104, 3623-3640.  This function was written, initially in MATLAB, in 2010. 
 33   
 34   
 35  Code origin 
 36  =========== 
 37   
 38  This is the model of the numerical solution for the 2-site Bloch-McConnell equations.  It originates as optimization function number 1 from the fitting_main_kex.py script from Mathilde Lescanne, Paul Schanda, and Dominique Marion (see U{http://thread.gmane.org/gmane.science.nmr.relax.devel/4138}, U{https://web.archive.org/web/https://gna.org/task/?7712#comment2} and U{https://web.archive.org/web/https://gna.org/support/download.php?file_id=18262}). 
 39   
 40   
 41  Links 
 42  ===== 
 43   
 44  More information on the NS CPMG 2-site 3D model can be found in the: 
 45   
 46      - U{relax wiki<http://wiki.nmr-relax.com/NS_CPMG_2-site_3D>}, 
 47      - U{relax manual<http://www.nmr-relax.com/manual/The_reduced_NS_2_site_3D_CPMG_model.html>}, 
 48      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_CPMG_2-site_3D>}. 
 49   
 50  More information on the NS CPMG 2-site 3D full model can be found in the: 
 51   
 52      - U{relax wiki<http://wiki.nmr-relax.com/NS_CPMG_2-site_3D_full>}, 
 53      - U{relax manual<http://www.nmr-relax.com/manual/The_full_NS_2_site_3D_CPMG_model.html>}, 
 54      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_CPMG_2-site_3D_full>}. 
 55  """ 
 56   
 57  # Python module imports. 
 58  from numpy import array, dot, fabs, float64, einsum, isfinite, log, min, multiply, rollaxis, sum 
 59  from numpy.linalg import matrix_power 
 60  from numpy.ma import fix_invalid, masked_where 
 61   
 62  # relax module imports. 
 63  from lib.float import isNaN 
 64  from lib.dispersion.matrix_exponential import matrix_exponential 
 65   
 66  # Repetitive calculations (to speed up calculations). 
 67  m_r10a = array([ 
 68      [0,  0,  0,  0,  0,  0,  0], 
 69      [0,  0,  0,  0,  0,  0,  0], 
 70      [0,  0,  0,  0,  0,  0,  0], 
 71      [1,  0,  0, -1,  0,  0,  0], 
 72      [0,  0,  0,  0,  0,  0,  0], 
 73      [0,  0,  0,  0,  0,  0,  0], 
 74      [0,  0,  0,  0,  0,  0,  0]], float64) 
 75   
 76  m_pA = array([ 
 77      [0,  0,  0,  0,  0,  0,  0], 
 78      [0,  0,  0,  0,  0,  0,  0], 
 79      [0,  0,  0,  0,  0,  0,  0], 
 80      [2,  0,  0,  0,  0,  0,  0], 
 81      [0,  0,  0,  0,  0,  0,  0], 
 82      [0,  0,  0,  0,  0,  0,  0], 
 83      [0,  0,  0,  0,  0,  0,  0]], float64) 
 84   
 85  m_r10b = array([ 
 86      [0,  0,  0,  0,  0,  0,  0], 
 87      [0,  0,  0,  0,  0,  0,  0], 
 88      [0,  0,  0,  0,  0,  0,  0], 
 89      [0,  0,  0,  0,  0,  0,  0], 
 90      [0,  0,  0,  0,  0,  0,  0], 
 91      [0,  0,  0,  0,  0,  0,  0], 
 92      [1,  0,  0,  0,  0,  0, -1]], float64) 
 93   
 94  m_pB = array([ 
 95      [0,  0,  0,  0,  0,  0,  0], 
 96      [0,  0,  0,  0,  0,  0,  0], 
 97      [0,  0,  0,  0,  0,  0,  0], 
 98      [0,  0,  0,  0,  0,  0,  0], 
 99      [0,  0,  0,  0,  0,  0,  0], 
100      [0,  0,  0,  0,  0,  0,  0], 
101      [2,  0,  0,  0,  0,  0,  0]], float64) 
102   
103  m_r20a = array([ 
104      [0,  0,  0,  0,  0,  0,  0], 
105      [0, -1,  0,  0,  0,  0,  0], 
106      [0,  0, -1,  0,  0,  0,  0], 
107      [0,  0,  0,  0,  0,  0,  0], 
108      [0,  0,  0,  0,  0,  0,  0], 
109      [0,  0,  0,  0,  0,  0,  0], 
110      [0,  0,  0,  0,  0,  0,  0]], float64) 
111   
112  m_r20b = array([ 
113      [0,  0,  0,  0,  0,  0,  0], 
114      [0,  0,  0,  0,  0,  0,  0], 
115      [0,  0,  0,  0,  0,  0,  0], 
116      [0,  0,  0,  0,  0,  0,  0], 
117      [0,  0,  0,  0, -1,  0,  0], 
118      [0,  0,  0,  0,  0, -1,  0], 
119      [0,  0,  0,  0,  0,  0,  0]], float64) 
120   
121  m_k_AB = array([ 
122      [0,  0,  0,  0,  0,  0,  0], 
123      [0, -1,  0,  0,  0,  0,  0], 
124      [0,  0, -1,  0,  0,  0,  0], 
125      [0,  0,  0, -1,  0,  0,  0], 
126      [0,  1,  0,  0,  0,  0,  0], 
127      [0,  0,  1,  0,  0,  0,  0], 
128      [0,  0,  0,  1,  0,  0,  0]], float64) 
129   
130  m_k_BA = array([ 
131      [0,  0,  0,  0,  0,  0,  0], 
132      [0,  0,  0,  0,  1,  0,  0], 
133      [0,  0,  0,  0,  0,  1,  0], 
134      [0,  0,  0,  0,  0,  0,  1], 
135      [0,  0,  0,  0, -1,  0,  0], 
136      [0,  0,  0,  0,  0, -1,  0], 
137      [0,  0,  0,  0,  0,  0, -1]], float64) 
138   
139  m_wA = array([ 
140      [0,  0,  0,  0,  0,  0,  0], 
141      [0,  0, -1,  0,  0,  0,  0], 
142      [0,  1,  0,  0,  0,  0,  0], 
143      [0,  0,  0,  0,  0,  0,  0], 
144      [0,  0,  0,  0,  0,  0,  0], 
145      [0,  0,  0,  0,  0,  0,  0], 
146      [0,  0,  0,  0,  0,  0,  0]], float64) 
147   
148  m_wB = array([ 
149      [0,  0,  0,  0,  0,  0,  0], 
150      [0,  0,  0,  0,  0,  0,  0], 
151      [0,  0,  0,  0,  0,  0,  0], 
152      [0,  0,  0,  0,  0,  0,  0], 
153      [0,  0,  0,  0,  0, -1,  0], 
154      [0,  0,  0,  0,  1,  0,  0], 
155      [0,  0,  0,  0,  0,  0,  0]], float64) 
156   
157   
158 -def rcpmg_3d_rankN(R1A=None, R1B=None, R2A=None, R2B=None, pA=None, pB=None, dw=None, k_AB=None, k_BA=None, tcp=None):
159 """Definition of the 3D exchange matrix, for rank [NE][NS][NM][NO][ND][7][7]. 160 161 @keyword R1A: The longitudinal, spin-lattice relaxation rate for state A. 162 @type R1A: numpy float array of rank [NE][NS][NM][NO][ND] 163 @keyword R1B: The longitudinal, spin-lattice relaxation rate for state B. 164 @type R1B: numpy float array of rank [NE][NS][NM][NO][ND] 165 @keyword R2A: The transverse, spin-spin relaxation rate for state A. 166 @type R2A: numpy float array of rank [NE][NS][NM][NO][ND] 167 @keyword R2B: The transverse, spin-spin relaxation rate for state B. 168 @type R2B: numpy float array of rank [NE][NS][NM][NO][ND] 169 @keyword pA: The population of state A. 170 @type pA: float 171 @keyword pB: The population of state B. 172 @type pB: float 173 @keyword dw: The chemical exchange difference between states A and B in rad/s. 174 @type dw: numpy float array of rank [NE][NS][NM][NO][ND] 175 @keyword k_AB: The forward exchange rate from state A to state B. 176 @type k_AB: float 177 @keyword k_BA: The reverse exchange rate from state B to state A. 178 @type k_BA: float 179 @keyword tcp: The tau_CPMG times (1 / 4.nu1). 180 @type tcp: numpy float array of rank [NE][NS][NM][NO][ND] 181 @return: The relaxation matrix. 182 @rtype: numpy float array of rank [NE][NS][NM][NO][ND][7][7] 183 """ 184 185 # The omega frequencies for states A and B (state A is assumed to be at zero frequency). 186 wA = 0.0 187 wB = dw 188 189 r10a_tcp = R1A * tcp 190 r10b_tcp = R1B * tcp 191 r20a_tcp = R2A * tcp 192 r20b_tcp = R2B * tcp 193 pA_tcp = pA * tcp 194 pB_tcp = pB * tcp 195 dw_tcp = dw * tcp 196 k_AB_tcp = k_AB * tcp 197 k_BA_tcp = k_BA * tcp 198 wA_tcp = wA * tcp 199 wB_tcp = wB * tcp 200 201 # Create the matrix. 202 # Multiply and expand. 203 m_r10a_tcp = multiply.outer( r10a_tcp, m_r10a ) 204 m_pA_tcp = multiply.outer( pA_tcp, m_pA ) 205 206 m_r10b_tcp = multiply.outer( r10b_tcp, m_r10b ) 207 m_pB_tcp = multiply.outer( pB_tcp, m_pB ) 208 209 m_r20a_tcp = multiply.outer( r20a_tcp, m_r20a ) 210 m_r20b_tcp = multiply.outer( r20b_tcp, m_r20b ) 211 212 m_k_AB_tcp = multiply.outer( k_AB_tcp, m_k_AB ) 213 m_k_BA_tcp = multiply.outer( k_BA_tcp, m_k_BA ) 214 215 m_wA_tcp = multiply.outer( wA_tcp, m_wA ) 216 m_wB_tcp = multiply.outer( wB_tcp, m_wB ) 217 218 # Collect matrix. 219 c_mat = (m_r10a_tcp * m_pA_tcp + m_r10b_tcp * m_pB_tcp 220 + m_r20a_tcp + m_r20b_tcp 221 + m_k_AB_tcp + m_k_BA_tcp 222 + m_wA_tcp + m_wB_tcp ) 223 224 # Return the matrix. 225 return c_mat
226 227
228 -def r2eff_ns_cpmg_2site_3D(r180x=None, M0=None, M0_T=None, r10a=0.0, r10b=0.0, r20a=None, r20b=None, pA=None, dw=None, dw_orig=None, kex=None, inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None):
229 """The 2-site numerical solution to the Bloch-McConnell equation. 230 231 This function calculates and stores the R2eff values. 232 233 234 @keyword r180x: The X-axis pi-pulse propagator. 235 @type r180x: numpy float64, rank-2, 7D array 236 @keyword M0: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 237 @type M0: numpy float array of rank [NE][NS][NM][NO][ND][7][1] 238 @keyword M0_T: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations, where the outer two axis has been swapped for efficient dot operations. 239 @type M0_T: numpy float array of rank [NE][NS][NM][NO][ND][1][7] 240 @keyword r10a: The R1 value for state A. 241 @type r10a: float 242 @keyword r10b: The R1 value for state B. 243 @type r10b: float 244 @keyword r20a: The R2 value for state A in the absence of exchange. 245 @type r20a: numpy float array of rank [NE][NS][NM][NO][ND] 246 @keyword r20b: The R2 value for state B in the absence of exchange. 247 @type r20b: numpy float array of rank [NE][NS][NM][NO][ND] 248 @keyword pA: The population of state A. 249 @type pA: float 250 @keyword dw: The chemical exchange difference between states A and B in rad/s. 251 @type dw: numpy float array of rank [NE][NS][NM][NO][ND] 252 @keyword dw_orig: The chemical exchange difference between states A and B in ppm. This is only for faster checking of zero value, which result in no exchange. 253 @type dw_orig: numpy float array of rank-1 254 @keyword kex: The kex parameter value (the exchange rate in rad/s). 255 @type kex: float 256 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds). 257 @type inv_tcpmg: numpy float array of rank [NE][NS][NM][NO][ND] 258 @keyword tcp: The tau_CPMG times (1 / 4.nu1). 259 @type tcp: numpy float array of rank [NE][NS][NM][NO][ND] 260 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies. 261 @type back_calc: numpy float array of rank [NE][NS][NM][NO][ND] 262 @keyword num_points: The number of points on the dispersion curve, equal to the length of the tcp and back_calc arguments. 263 @type num_points: numpy int array of rank [NE][NS][NM][NO] 264 @keyword power: The matrix exponential power array. 265 @type power: numpy int array of rank [NE][NS][NM][NO][ND] 266 """ 267 268 # Flag to tell if values should be replaced if math function is violated. 269 t_dw_zero = False 270 271 # Catch parameter values that will result in no exchange, returning flat R2eff = R20 lines (when kex = 0.0, k_AB = 0.0). 272 if pA == 1.0 or kex == 0.0: 273 back_calc[:] = r20a 274 return 275 276 # Test if dw is zero. Create a mask for the affected spins to replace these with R20 at the end of the calculationWait for replacement, since this is spin specific. 277 if min(fabs(dw_orig)) == 0.0: 278 t_dw_zero = True 279 mask_dw_zero = masked_where(dw == 0.0, dw) 280 281 # Once off parameter conversions. 282 pB = 1.0 - pA 283 k_BA = pA * kex 284 k_AB = pB * kex 285 286 # This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 287 M0_T[:, :, :, :, :, 0, 1] = pA 288 M0_T[:, :, :, :, :, 0, 4] = pB 289 M0[:, :, :, :, :, 1, 0] = pA 290 M0[:, :, :, :, :, 4, 0] = pB 291 292 # Extract the total numbers of experiments, number of spins, number of magnetic field strength, number of offsets, maximum number of dispersion point. 293 NE, NS, NM, NO, ND = back_calc.shape 294 295 # The matrix R that contains all the contributions to the evolution, i.e. relaxation, exchange and chemical shift evolution. 296 R_mat = rcpmg_3d_rankN(R1A=r10a, R1B=r10b, R2A=r20a, R2B=r20b, pA=pA, pB=pB, dw=dw, k_AB=k_AB, k_BA=k_BA, tcp=tcp) 297 298 # This matrix is a propagator that will evolve the magnetization with the matrix R for a delay tcp. 299 Rexpo_mat = matrix_exponential(R_mat) 300 301 # The the essential evolution matrix. 302 # This is a dot product of the outer [7][7] matrix of the Rexpo_mat and r180x matrixes, which 303 # have the shape [NE][NS][NM][NO][ND][7][7] and [7][7]. 304 # This can be achieved by using numpy einsum, and where ellipsis notation will use the last axis. 305 evolution_matrix_mat = einsum('...ij,...jk', Rexpo_mat, r180x) 306 evolution_matrix_mat = einsum('...ij,...jk', evolution_matrix_mat, Rexpo_mat) 307 evolution_matrix_mat = einsum('...ij,...jk', evolution_matrix_mat, evolution_matrix_mat) 308 309 # Roll axis around. 310 evolution_matrix_T_mat = rollaxis(evolution_matrix_mat, 6, 5) 311 312 # Preform the initial magnetisation. 313 evolution_matrix_T_M0_mat = einsum('...ij,...jk', M0_T, evolution_matrix_T_mat) 314 315 # Loop over the spins 316 for si in range(NS): 317 # Loop over the spectrometer frequencies. 318 for mi in range(NM): 319 # Extract number of points. 320 num_points_si_mi = int(num_points[0, si, mi, 0]) 321 322 # Loop over the time points, back calculating the R2eff values. 323 for di in range(num_points_si_mi): 324 # Extract the values from the higher dimensional arrays. 325 inv_tcpmg_si_mi_di = inv_tcpmg[0, si, mi, 0, di] 326 power_si_mi_di = int(power[0, si, mi, 0, di]) 327 r20a_si_mi_di = r20a[0, si, mi, 0, di] 328 329 # Initial magnetisation. 330 Mint_T_i = evolution_matrix_T_M0_mat[0, si, mi, 0, di] 331 332 # This matrix is a propagator that will evolve the magnetization with the matrix R for a delay tcp. 333 evolution_matrix_T_i = evolution_matrix_T_mat[0, si, mi, 0, di] 334 335 # Get which power to raise the matrix to. 336 l = int(power_si_mi_di-1) 337 338 # Raise the square evolution matrix to the power l. 339 evolution_matrix_T_power_i = matrix_power(evolution_matrix_T_i, l) 340 341 # Evolve the magnetisation. 342 Mint_T_i = dot(Mint_T_i, evolution_matrix_T_power_i) 343 344 # The next lines calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential. 345 Mx = Mint_T_i[0][1] / pA 346 if Mx <= 0.0 or isNaN(Mx): 347 back_calc[0, si, mi, 0, di] = r20a_si_mi_di 348 else: 349 back_calc[0, si, mi, 0, di] = - inv_tcpmg_si_mi_di * log(Mx) 350 351 # Replace data in array. 352 # If dw is zero. 353 if t_dw_zero: 354 back_calc[mask_dw_zero.mask] = r20a[mask_dw_zero.mask] 355 356 # Catch errors, taking a sum over array is the fastest way to check for 357 # +/- inf (infinity) and nan (not a number). 358 if not isfinite(sum(back_calc)): 359 # Replaces nan, inf, etc. with fill value. 360 fix_invalid(back_calc, copy=False, fill_value=1e100)
361