Hi,
fg is the frequency of the ground state (call it "A" or "G"), while fe is the excited state frequency.
Is this the chemical shift in Hz? Or the chemical shift difference in Hz? The code does not use the 'fe' parameter, so how is the chemical shift of the excited state handled? Does this mean that chemical shifts need to be explicitly given? How does this work with intermediate to fast exchange where the shift of the peak is a mixture of the two chemical shifts? I'm still not sure how to handle this in relax. I may need to add a new dispersion parameter(s) to the definitions at the start of the Relax_disp class in specific_analyses/relax_disp/api.py. Though very easy to expand, this currently only supports the following parameters for chemical shifts: phi_ex = pA.pB.dw**2 (in ppm^2) padw2 = pA.dw**2 (in ppm^2) dw = (in ppm) where dw is the chemical shift difference. In all of the analytic models, only dw is used.
tcpmg is the total duration of the CPMG element. At some instances this is called Tc instead of tcpmg. Not to be confused with the delay before and after each pulse, which is tcp. I.e., between two 180deg pulses there are two delays tcp.
Thank you for the additional details. I have now updated the lib/dispersion/ns_2site_star.py file with the details (http://article.gmane.org/gmane.science.nmr.relax.scm/18046).
I saw that you keep the convention of "G" and "E" in the functions. Is this also used in the existing functions, or do you use "A" and "B" there? It might be cleaner to have a coherent naming convention.
I am in the process of converting. I think I'll use the arbitrary mapping: G -> A E -> B This is for consistency within relax between this code and the analytic models, as well as to minimise user confusion. It will also allow relax to more logically support a C state for the 3-site models.
Here are some comments to the function: 70 # Initialise some structures. 71 Rr = -1.0 * matrix([[r20b, 0.0],[0.0, r20a]]) 72 Rex = -1.0 * matrix([[kge, -keg],[-kge, keg]]) 73 RCS = complex(0.0, -1.0) * matrix([[0.0, 0.0],[0.0, fg]]) 74 R = Rr + Rex + RCS 75 cR2 = conj(R) * 2.0 76 77 kex = kge + keg # conversion of kinetic rates 78 IGeq = keg / kex # calculate relative populations - will be used for generating the equilibrium magnetizations 79 IEeq = kge / kex # like preceding line, but for excited state 80 M0 = matrix([[IGeq], [IEeq]]) # This is a vector that contains the initial magnetizations, corresponding to ground (G) and excited (E) state transverse magnetizations 81 82 # Replicated calculations for faster operation. 83 inv_tcpmg = 1.0 / tcpmg 84 85 # Loop over the time points, back calculating the R2eff values. 86 for i in range(num_points): 87 # Replicated calculations for faster operation. 88 tcp = 0.25 / cpmg_frqs[i] 89 expm_R_tcp = expm(R*tcp) # This matrix is a propagator that will evolve the magnetization with the matrix R for a delay tcp 90 91 prop_2 = dot(dot(expm_R_tcp, expm(cR2*tcp)), expm_R_tcp) # This is the propagator for an element of delay tcp - 180degpulse -2 times delay tcp - 180degpulse - delay tau, i.e. for 2 times tau-180-tau 92 93 prop_total = square_matrix_power(prop_2, cpmg_frqs[i]*tcpmg) # Now create the total propagator, that will evolve the magnetization under the CPMG train, ie. it applies the above tau-180-tau-tau-180-tau so many times as required for the CPMG frequency under consideration. 94 95 Moft = prop_total * M0 # 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 96 97 Mgx = Moft[0].real / M0[0] # This and the next line calculate the R2eff, using a two-point approximation, i.e. assuming that the decay is mono-exponential. 98 back_calc[i]= -inv_tcpmg * log(Mgx)
These comments are awesome!!! It's far better than what I could have added. The code in lib/dispersion/ns_2site_star.py is now a perfect teaching aid that can be shown directly to students. It will also allow NMR spectroscopists who are not experts to quickly learn from the code. I have added all of these comments and updated your copyright notice appropriately.
Thanks for the implementation!
Thank you for contributing the code. I hope this will up and running soon and, together with the tools built into relax, that the calculations will even be faster. Cheers, Edward