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 12, 2013 - 09:43:
Thanks!  I've updated the code with this information.  See commit
r20276 (http://article.gmane.org/gmane.science.nmr.relax.scm/18032 -
this link takes some time to be generated) and the code at
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/ns_2site_star.py?revision=20276&view=markup
after clicking on 'as text' (or type 'svn up' in your checked out copy
of the relax_disp branch).

In relax, the pA parameter together with kex is used.  This is the
current convention and the reason is quite important - it due to how
relax performs optimisation for the relaxation dispersion analysis.
For optimisation, the simplex algorithm is used
(http://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method) together
with the log-barrier constraint algorithm
(http://en.wikipedia.org/wiki/Barrier_function#Logarithmic_barrier_function).
 This allows for the fastest and most accurate optimisation possible
when gradients and Hessians cannot be derived (1st and 2nd partial
derivatives), as is the case for the numerical solutions.  I hope we
can get this model up and running soon to see how fast we can make it
run.

The important part about having pA as a parameter is that I can use
the constraints 0.5 <= pA <= 1.0.  The 0.5 limit cuts the optimisation
space in half.  Both sides of this limit are mirror image spaces, so
it is numerically inefficient to search both spaces.  The 1.0 limit is
also important as numerically pA will often go above 1 during
optimisation.  I have seen this in a number of cases.  Sometimes pA
comes back, but in most cases it runs off into impossibly high values
(compensated by negative pB values).  The result is incredibly slow
optimisation to a nonsensical solution for that spin cluster.  This
occurs when the model is not suitable for the data.  Without a pA
parameter, I cannot stop this behaviour.

If you could advise how to change from {kge, keg} to {pA, kex} in your
code, it would be appreciated.

Cheers,

Edward


On 12 July 2013 07:51, Paul Schanda <paul.schanda@xxxxxx> wrote:
Hi Edward,

Yes, R2E and R2G are the exchange-free R2 values. Sometimes also referred to
as R20 or R2infinity. They are the Redfield-relaxation part.
Currently, the functions do not use populations (pA, pB) but they use the
two exchange rates (forward, backward), rather than one rate constant and
one population. The two different approaches are in essence the same, as
populations can be calculated from forward/backward rates in a
straightforward manner.
We can easily change this if it is more in the logic of the existing code.

paul


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
paul.schanda@xxxxxx





_______________________________________________
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




Related Messages


Powered by MHonArc, Updated Fri Jul 12 11:40:07 2013