Hi Edward,
fg is the frequency of the ground state (call it "A" or "G"), while fe is the excited state frequency. 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.
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.
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) |
Thanks for the implementation!
Paul
Hi Paul and Mathilde,
I have attached this file to the task report #7712 (https://gna.org/task/?7712) as a permanent reference (https://gna.org/support/download.php?file_id=18263). For new files it would be appreciated if they were attached to this same task with comments rather than sent to the public mailing list, as this avoids unnecessary strain on the open source infrastructure, and does not require a Gna! login. Cheers.
I have now added the parameter conversions in that file to the func_ns_2site_star() function of the target_functions/relax_disp.py file, instead of in the r2eff_ns_2site_star() function of lib/dispersion/ns_2site_star.py. This is to minimise the mathematical operations for faster code. You can see it in revision r20284. I may do this for kex, IGeq, IEeq and M0 as well once the code is functional. This will cause large speed ups if large clusters of spins are used or data at many magnetic field strengths are collected.
I still have two problems - the 'fg' and 'tcpmg' parameters. Could you explain these? I'm pretty sure I know what they are, but it would be good to confirm. Another thing that would be nice to add to the code in lib/dispersion/ns_2site_star.py would be a short description of the matrices and what the operations are. I.e. a simple comment for what each data structure is.
I have seen that the mpower() function is in this file, so I'll just copy that directly somewhere into the relax library.
Regards,
Edward
On 12 July 2013 12:09, Paul Schanda <paul.schanda@xxxxxx> wrote:
Hi Edward,
I understand your point of fitting pA and kex instead of keg and kge (or kAB/kBA). Please find attached a program with the convention of kex and pB. The states A and B are called G and E (ground/excited), so the parameters are generally associated with E and G instead of A and B.
You will find in the attached program six different implementations, that can be chosen with the variable a set to 1-6. In practice, however, we are using only the "2D complex" numerical approach (a=6), or the Maple-derived approach (a=5). The latter is faster, and is most-often used. It has not been published, though. The numerical Bloch-McConnell treatment is described in many introductory NMR books, and also in e.g. McConnell, H. J Chem Phys 1958, 28, 430–431. (the original reference) Palmer, A. G. Chem Rev 2004, 104, 3623–3640. Palmer, A. G.; Massi, F. Chem Rev 2006, 106, 1700–1719. Baldwin, A. J.; Kay, L. E. J Biomol NMR 2013. The latter two are R1rho, but the formalism of Bloch-McConnell is the same.
Paul
Hi Paul,
I have now implemented step 2 of the tutorial for adding relaxation dispersion models (http://article.gmane.org/gmane.science.nmr.relax.devel/3907, in the commit r20277 at http://article.gmane.org/gmane.science.nmr.relax.scm/18033). This is one aspect that the relax user sees. I have assumed that 'dw' is a parameter of this model - i.e. the Delta_omega chemical shift difference between states A and B. But I am not sure if this is used in your code and if it relates to the fg parameter. Would you know the relationship?
Cheers,
Edward
On 11 July 2013 18:23, Edward d'Auvergne <edward@xxxxxxxxxxxxx> wrote:
Hello,
Another question I have is what are the parameters of this model? It relates to the function arguments currently labelled as unknown in the function docstring. I would like to implement step 2 of the tutorial on adding relaxation dispersion models to relax (the relax_disp.select_model user function front end, see http://article.gmane.org/gmane.science.nmr.relax.devel/3907), but I'm not sure how to list the parameters of the model for the user. Is there a population parameter (pA or pB)? Are R2E and R2G the same as the R20A and R20B rates (the R2 rates in the absence of exchange for states A and B for a fixed magnetic field strength)?
Cheers,
Edward
On 11 July 2013 16:14, Edward d'Auvergne <edward@xxxxxxxxxxxxx> wrote:
Hi Paul,
I have now incorporated this code into relax (see http://article.gmane.org/gmane.science.nmr.relax.scm/18025). Could you please check the code, comments, and docstrings carefully and see if there are any issues? Could you also tell me if anyone else should be included in the copyright notice? It would be useful to replace the 'Unknown' elements in the function docstring with basic descriptions of the parameter and its Python object type. And any suggestions for comments describing in basic NMR terms what the code is doing (the physics of the propagations, for example) would also be appreciated.
You can get the new code by typing the following in the base directory of the checked out copy of the relax_disp branch:
$ svn up
There was an indentation problem that I have tried to solve. Also I cannot determine where the mpower() function comes from - it appears not to be part of numpy or scipy. Is this an outside function you wrote? I have also made a number of significant speed ups of the code (see article.gmane.org/gmane.science.nmr.relax.scm/18029). Note that this cannot be selected as a model in the dispersion analysis yet.
Cheers,
Edward
On 11 July 2013 13:47, Paul Schanda <paul.schanda@xxxxxx> wrote:
Hi Edward,
Here is a function that does a numerical fit of CPMG. It 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. It returns calculated R2eff values, so it can be used for numerical fitting of CPMG. It is an addition to the functions that I previously sent to you. The approach of Bloch-McConnell can be found in chapter 3.1 of Palmer, A. G. Chem Rev 2004, 104, 3623–3640. I wrote this function (initially in MATLAB) in 2010. I agree that this code is released under the GPL3 or higher licence.
Paul
def func1(R2E,R2G,fg,kge,keg, cpmgfreq,tcpmg): kex=kge+keg Rcalc_array=np.zeros(len(cpmgfreq)) Rr = -1.*np.matrix([[R2G, 0.],[0., R2E]]) Rex = -1.*np.matrix([[kge,-keg],[-kge, keg]]) RCS = complex(0.,-1.)*np.matrix([[0. ,0.],[0.,fg]]) R = Rr + Rex + RCS cR = np.conj(R)
IGeq=keg/kex; IEeq=kge/kex M0=np.matrix([[IGeq],[IEeq]])
for k in range(len(cpmgfreq)): tcp=1./(4.*cpmgfreq[k]) prop_2 =np.dot(np.dot(sci.expm(R*tcp),sci.expm(cR*2.*tcp)),sci.expm(R*tcp)) prop_total =mpower(prop_2,cpmgfreq[k]*tcpmg)
Moft = prop_total*M0
Mgx=Moft[0].real/M0[0] Rcalc=-(1./tcpmg)*math.log(Mgx); Rcalc_array[k]=Rcalc return Rcalc_array
_______________________________________________ relax (http://www.nmr-relax.com)
This is the relax-devel mailing list relax-devel@xxxxxxx
To unsubscribe from this list, get a password reminder, or change your subscription options, visit the list information page at https://mail.gna.org/listinfo/relax-devel
***************************
Paul Schanda
|