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