mailRe: numerical cpmg fit


Others Months | Index by Date | Thread Index
>>   [Date Prev] [Date Next] [Thread Prev] [Thread Next]

Header


Content

Posted by Edward d'Auvergne on July 17, 2013 - 11:11:
Cheers, I have now updated all the code to match and to have comments
about this.

Regards,

Edward


On 16 July 2013 08:01, Paul Schanda <paul.schanda@xxxxxx> wrote:
Yes, that's right. As one of the states is assumed to be at zero frequency,
the frequency of the other state is the same as the chemical shift
difference.

paul


Am 15.07.2013 um 10:56 schrieb "Edward d'Auvergne" <edward@xxxxxxxxxxxxx>:

Hi,

Could it be that the 'fg' parameter in the original code at
http://thread.gmane.org/gmane.science.nmr.relax.devel/4132 is actually
the chemical shift difference?  Is 'fe' assumed to be zero, hence 'fg
= fe + df', where 'df' is the chemical shift difference in rad/second?

Regards,

Edward


On 15 July 2013 10:01, Edward d'Auvergne <edward@xxxxxxxxxxxxx> wrote:

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




***************************

Paul Schanda
paul.schanda@xxxxxx







Related Messages


Powered by MHonArc, Updated Wed Jul 17 11:40:06 2013