lib.dispersion.ns_cpmg_2site

1 ############################################################################### 2 # # 3 # Copyright (C) 2000-2001 Nikolai Skrynnikov # 4 # Copyright (C) 2000-2001 Martin Tollinger # 5 # Copyright (C) 2010-2013 Paul Schanda # 6 # Copyright (C) 2013 Mathilde Lescanne # 7 # Copyright (C) 2013 Dominique Marion # 8 # Copyright (C) 2013-2014 Edward d'Auvergne # 9 # Copyright (C) 2014 Troels E. Linnet # 10 # # 11 # This file is part of the program relax (http://www.nmr-relax.com). # 12 # # 13 # This program is free software: you can redistribute it and/or modify # 14 # it under the terms of the GNU General Public License as published by # 15 # the Free Software Foundation, either version 3 of the License, or # 16 # (at your option) any later version. # 17 # # 18 # This program is distributed in the hope that it will be useful, # 19 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 20 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 21 # GNU General Public License for more details. # 22 # # 23 # You should have received a copy of the GNU General Public License # 24 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 25 # # 26 ############################################################################### 27 28 # Module docstring. 29 """The numerical fit of 2-site Bloch-McConnell equations for CPMG-type experiments, the U{NS CPMG 2-site star<http://wiki.nmr-relax.com/NS_CPMG_2-site_star>} and U{NS CPMG 2-site star full<http://wiki.nmr-relax.com/NS_CPMG_2-site_star_full>} models. 30 31 Description 32 =========== 33 34 The function uses an explicit matrix that contains relaxation, exchange and chemical shift terms. It does the 180deg pulses in the CPMG train with conjugate complex matrices. The approach of Bloch-McConnell can be found in chapter 3.1 of Palmer, A. G. 2004 I{Chem. Rev.}, B{104}, 3623-3640. This function was written, initially in MATLAB, in 2010. 35 36 37 Code origin 38 =========== 39 40 The code was submitted at U{http://thread.gmane.org/gmane.science.nmr.relax.devel/4132} by Paul Schanda. 41 42 43 Links 44 ===== 45 46 More information on the NS CPMG 2-site star model can be found in the: 47 48 - U{relax wiki<http://wiki.nmr-relax.com/NS_CPMG_2-site_star>}, 49 - U{relax manual<http://www.nmr-relax.com/manual/The_reduced_NS_2_site_star_CPMG_model.html>}, 50 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_CPMG_2-site_star>}. 51 52 More information on the NS CPMG 2-site star full model can be found in the: 53 54 - U{relax wiki<http://wiki.nmr-relax.com/NS_CPMG_2-site_star_full>}, 55 - U{relax manual<http://www.nmr-relax.com/manual/The_full_NS_2_site_star_CPMG_model.html>}, 56 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_CPMG_2-site_star_full>}. 57 """ 58 59 # Python module imports. 60 from numpy import add, array, conj, dot, einsum, fabs, float64, isfinite, log, min, multiply, sum 61 from numpy.ma import fix_invalid, masked_where 62 from numpy.linalg import matrix_power 63 64 # relax module imports. 65 from lib.float import isNaN 66 from lib.dispersion.matrix_exponential import matrix_exponential 67 68 # Repetitive calculations (to speed up calculations). 69 m_r20a = array([ 70 [-1, 0], 71 [ 0, 0]], float64) 72 73 m_r20b = array([ 74 [ 0, 0], 75 [ 0, -1]], float64) 76 77 m_k_AB = array([ 78 [-1, 0], 79 [ 1, 0]], float64) 80 81 m_k_BA = array([ 82 [ 0, 1], 83 [ 0, -1]], float64) 84 85 m_dw = array([ 86 [ 0, 0], 87 [ 0, 1]], float64) 88 89

90 -def rcpmg_star_rankN(R2A=None, R2B=None, dw=None, k_AB=None, k_BA=None, tcp=None):

91 """Definition of the exchange matrix, for rank [NE][NS][NM][NO][ND][2][2]. 92 93 @keyword R2A: The transverse, spin-spin relaxation rate for state A. 94 @type R2A: numpy float array of rank [NE][NS][NM][NO][ND] 95 @keyword R2B: The transverse, spin-spin relaxation rate for state B. 96 @type R2B: numpy float array of rank [NE][NS][NM][NO][ND] 97 @keyword dw: The chemical exchange difference between states A and B in rad/s. 98 @type dw: numpy float array of rank [NE][NS][NM][NO][ND] 99 @keyword k_AB: The forward exchange rate from state A to state B. 100 @type k_AB: float 101 @keyword k_BA: The reverse exchange rate from state B to state A. 102 @type k_BA: float 103 @keyword tcp: The tau_CPMG times (1 / 4.nu1). 104 @type tcp: numpy float array of rank [NE][NS][NM][NO][ND] 105 @return: The relaxation matrix R and complex conjugate cR2. 106 @rtype: numpy float array of rank [NE][NS][NM][NO][ND][2][2] 107 """ 108 109 # Pre-multiply with tcp. 110 r20a_tcp = R2A * tcp 111 r20b_tcp = R2B * tcp 112 k_AB_tcp = k_AB * tcp 113 k_BA_tcp = k_BA * tcp 114 # Complex dw. 115 dw_tcp_C = dw * tcp * -1j 116 117 # Create matrix for collection of Rr matrix. 118 # The matrix that contains only the R2 relaxation terms ("Redfield relaxation", i.e. non-exchange broadening). 119 #Rr[0, 0] = -R2A_si_mi 120 #Rr[1, 1] = -R2B_si_mi 121 122 # Multiply and expand. 123 m_r20a_tcp = multiply.outer( r20a_tcp, m_r20a ) 124 m_r20b_tcp = multiply.outer( r20b_tcp, m_r20b ) 125 126 # Collect Rr matrix. 127 Rr_mat = (m_r20a_tcp + m_r20b_tcp) 128 129 # Create matrix for collection of Rex. 130 # Set up the matrix that contains the exchange terms between the two states A and B. 131 #Rex[0, 0] = -k_AB 132 #Rex[0, 1] = k_BA 133 #Rex[1, 0] = k_AB 134 #Rex[1, 1] = -k_BA 135 136 # Multiply and expand. 137 m_k_AB_tcp = multiply.outer( k_AB_tcp, m_k_AB ) 138 m_k_BA_tcp = multiply.outer( k_BA_tcp, m_k_BA ) 139 140 # Collect Rex matrix. 141 Rex_mat = (m_k_AB_tcp + m_k_BA_tcp) 142 143 # Create the matrix for RCS. 144 # The matrix that contains the chemical shift evolution. It works here only with X magnetization, and the complex notation allows to evolve in the transverse plane (x, y). The chemical shift for state A is assumed to be zero. 145 #RCS[1, 1] = complex(0.0, -dw_si_mi) 146 147 # Multiply and expand. 148 m_dw_tcp_C = multiply.outer( dw_tcp_C, m_dw ) 149 150 # Collect RCS matrix. 151 RCS_mat = m_dw_tcp_C 152 153 # The matrix R that contains all the contributions to the evolution, i.e. relaxation, exchange and chemical shift evolution. 154 R_mat = add(Rr_mat, Rex_mat) 155 R_mat = add(R_mat, RCS_mat) 156 157 # This is the complex conjugate of the above. It allows the chemical shift to run in the other direction, i.e. it is used to evolve the shift after a 180 deg pulse. The factor of 2 is to minimise the number of multiplications for the prop_2 matrix calculation. 158 cR2_mat = conj(R_mat) * 2.0 159 160 # Return the matrixes. 161 return R_mat, cR2_mat, Rr_mat, Rex_mat, RCS_mat

162 163

164 -def r2eff_ns_cpmg_2site_star(M0=None, 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):

165 """The 2-site numerical solution to the Bloch-McConnell equation using complex conjugate matrices. 166 167 This function calculates and stores the R2eff values. 168 169 170 @keyword M0: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 171 @type M0: numpy float64, rank-1, 2D array 172 @keyword r20a: The R2 value for state A in the absence of exchange. 173 @type r20a: numpy float array of rank [NE][NS][NM][NO][ND] 174 @keyword r20b: The R2 value for state B in the absence of exchange. 175 @type r20b: numpy float array of rank [NE][NS][NM][NO][ND] 176 @keyword pA: The population of state A. 177 @type pA: float 178 @keyword dw: The chemical exchange difference between states A and B in rad/s. 179 @type dw: numpy float array of rank [NE][NS][NM][NO][ND] 180 @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. 181 @type dw_orig: numpy float array of rank-1 182 @keyword kex: The kex parameter value (the exchange rate in rad/s). 183 @type kex: float 184 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds). 185 @type inv_tcpmg: numpy float array of rank [NE][NS][NM][NO][ND] 186 @keyword tcp: The tau_CPMG times (1 / 4.nu1). 187 @type tcp: numpy float array of rank [NE][NS][NM][NO][ND] 188 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies. 189 @type back_calc: numpy float array of rank [NE][NS][NM][NO][ND] 190 @keyword num_points: The number of points on the dispersion curve, equal to the length of the tcp and back_calc arguments. 191 @type num_points: numpy int array of rank [NE][NS][NM][NO] 192 @keyword power: The matrix exponential power array. 193 @type power: numpy int array of rank [NE][NS][NM][NO][ND] 194 """ 195 196 # Flag to tell if values should be replaced if math function is violated. 197 t_dw_zero = False 198 199 # Catch parameter values that will result in no exchange, returning flat R2eff = R20 lines (when kex = 0.0, k_AB = 0.0). 200 if pA == 1.0 or kex == 0.0: 201 back_calc[:] = r20a 202 return 203 204 # 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. 205 if min(fabs(dw_orig)) == 0.0: 206 t_dw_zero = True 207 mask_dw_zero = masked_where(dw == 0.0, dw) 208 209 # Once off parameter conversions. 210 pB = 1.0 - pA 211 k_BA = pA * kex 212 k_AB = pB * kex 213 214 # This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 215 M0[0] = pA 216 M0[1] = pB 217 218 # Extract the total numbers of experiments, number of spins, number of magnetic field strength, number of offsets, maximum number of dispersion point. 219 NE, NS, NM, NO, ND = back_calc.shape 220 221 # The matrix R that contains all the contributions to the evolution, i.e. relaxation, exchange and chemical shift evolution. 222 R_mat, cR2_mat, Rr_mat, Rex_mat, RCS_mat = rcpmg_star_rankN(R2A=r20a, R2B=r20b, dw=dw, k_AB=k_AB, k_BA=k_BA, tcp=tcp) 223 224 # The the essential evolution matrix. 225 # This matrix is a propagator that will evolve the magnetization with the matrix R for a delay tcp. 226 eR_mat = matrix_exponential(R_mat) 227 ecR2_mat = matrix_exponential(cR2_mat) 228 229 # Preform the matrix. 230 # This is the propagator for an element of [delay tcp; 180 deg pulse; 2 times delay tcp; 180 deg pulse; delay tau], i.e. for 2 times tau-180-tau. 231 prop_2_mat = evolution_matrix_mat = einsum('...ij, ...jk', eR_mat, ecR2_mat) 232 prop_2_mat = evolution_matrix_mat = einsum('...ij, ...jk', prop_2_mat, eR_mat) 233 234 # Loop over the spins 235 for si in range(NS): 236 # Loop over the spectrometer frequencies. 237 for mi in range(NM): 238 # Extract the values from the higher dimensional arrays. 239 num_points_si_mi = int(num_points[0, si, mi, 0]) 240 241 # Loop over the time points, back calculating the R2eff values. 242 for di in range(num_points_si_mi): 243 # Extract the values from the higher dimensional arrays. 244 power_si_mi_di = int(power[0, si, mi, 0, di]) 245 246 # This is the propagator for an element of [delay tcp; 180 deg pulse; 2 times delay tcp; 180 deg pulse; delay tau], i.e. for 2 times tau-180-tau. 247 prop_2_i = prop_2_mat[0, si, mi, 0, di] 248 249 # Now create the total propagator that will evolve the magnetization under the CPMG train, i.e. it applies the above tau-180-tau-tau-180-tau so many times as required for the CPMG frequency under consideration. 250 prop_total = matrix_power(prop_2_i, power_si_mi_di) 251 252 # Now we apply the above propagator to the initial magnetization vector - resulting in the magnetization that remains after the full CPMG pulse train. It is called M of t (t is the time after the CPMG train). 253 Moft = dot(prop_total, M0) 254 255 # The next lines calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential. 256 Mx = Moft[0].real / M0[0] 257 if Mx <= 0.0 or isNaN(Mx): 258 back_calc[0, si, mi, 0, di] = 1e99 259 else: 260 back_calc[0, si, mi, 0, di]= -inv_tcpmg[0, si, mi, 0, di] * log(Mx) 261 262 # Replace data in array. 263 # If dw is zero. 264 if t_dw_zero: 265 back_calc[mask_dw_zero.mask] = r20a[mask_dw_zero.mask] 266 267 # Catch errors, taking a sum over array is the fastest way to check for 268 # +/- inf (infinity) and nan (not a number). 269 if not isfinite(sum(back_calc)): 270 # Replaces nan, inf, etc. with fill value. 271 fix_invalid(back_calc, copy=False, fill_value=1e100)

272

Source Code for Module lib.dispersion.ns_cpmg_2site_star