Hi Troels, If you need tips for hugely optimising this code, just ask ;) There is quite a lot that can be done. There are also a lot of changes which can be made to make the code fit better into relax - using the same parameter and variable naming conventions as in the other lib.dispersion modules, spacing around Python operators, spacing after commas (as the 2to3 program insists on), importing the numpy functions rather than using the module directly, etc.
Regards, Edward On 3 May 2014 21:35, <[email protected]> wrote: > Author: tlinnet > Date: Sat May 3 21:35:09 2014 > New Revision: 22940 > > URL: http://svn.gna.org/viewcvs/relax?rev=22940&view=rev > Log: > Implemented model B14 in the relax library. > > sr #3154: (https://gna.org/support/?3154) Implementation of Baldwin (2014) > B14 model - 2-site exact solution model for all time scales. > > This follows the tutorial for adding relaxation dispersion models at: > Tutorial_for_adding_relaxation_dispersion_models_to_relax#The_relax_library > > The code is raw implemented, with no optimisation. > > This is merely to test, that the spin parameters that created R2eff data, can > be found again after > grid searh and minimisation. > > The B14 model is explained in: http://wiki.nmr-relax.com/B14. > > Modified: > trunk/lib/dispersion/b14.py > > Modified: trunk/lib/dispersion/b14.py > URL: > http://svn.gna.org/viewcvs/relax/trunk/lib/dispersion/b14.py?rev=22940&r1=22939&r2=22940&view=diff > ============================================================================== > --- trunk/lib/dispersion/b14.py (original) > +++ trunk/lib/dispersion/b14.py Sat May 3 21:35:09 2014 > @@ -43,32 +43,32 @@ > The equation used is:: > > R2A0 + R2B0 + kex Ncyc 1 ( 1+y > 1-y ) > - R2eff = ------------------ - ------ * cosh^-1 * v1c - ------ ln( --- + > ------------------ * (v2 + 2*kAB*pD ) ) > + R2eff = ------------------ - ------ * cosh^-1 * v1c - ------ ln( --- + > ------------------ * (v2 + 2*kAB*pD ) ) , > 2 Trel Trel ( 2 > 2*sqrt(v1c^2 -1 ) ) > > 1 ( 1+y 1-y > ) > - = R2eff(CR72) - ------ ln( --- + ------------------ * (v2 + > 2*kAB*pD ) ) > + = R2eff(CR72) - ------ ln( --- + ------------------ * (v2 + > 2*kAB*pD ) ) , > Trel ( 2 2*sqrt(v1c^2 -1 ) > ) > > Which have these following definitions:: > > - v1c = F0 * cosh(tauCP * E0)- F2 * cosh(tauCP * E2) > - v1s = F0 * sinh(tauCP * E0)- F2 * sinh(tauCP * E2) > - v2*N = v1s * (OB-OA) + 4OB * F1^a * sinh(tauCP * E1) > - pD N = v1s + (F1^a + F1^b) * sinh(tauCP * E1) > - v3 = ( v2^2 + 4 * kBA * kAB * pD^2 )^1/2 > - y = ( (v1c-v3)/(v1c+v3) )^NCYC > + v1c = F0 * cosh(tauCP * E0)- F2 * cosh(tauCP * E2) , > + v1s = F0 * sinh(tauCP * E0)- F2 * sinh(tauCP * E2) , > + v2*N = v1s * (OB-OA) + 4OB * F1^a * sinh(tauCP * E1) , > + pD N = v1s + (F1^a + F1^b) * sinh(tauCP * E1) , > + v3 = ( v2^2 + 4 * kBA * kAB * pD^2 )^1/2 , > + y = ( (v1c-v3)/(v1c+v3) )^NCYC , > > Note, E2 is complex. If |x| denotes the complex modulus:: > > - cosh(tauCP * E2) = cos(tauCP * |E2|) > - sinh(tauCP * E2) = i sin(tauCP * |E2|) > + cosh(tauCP * E2) = cos(tauCP * |E2|) , > + sinh(tauCP * E2) = i sin(tauCP * |E2|) , > > The term pD is based on product of the off diagonal elements in the CPMG > propagator (Supplementary Section 3). > > It is interesting to consider the region of validity of the Carver Richards > result. The two results are equal when the correction is zero, which is true > when:: > > - sqrt(v1c^2-1) ~ v2 + 2*kAB * pD > + sqrt(v1c^2-1) ~ v2 + 2*kAB * pD , > > This occurs when 2*kAB * pD tends to zero, and so v2=v3. Setting "kAB * pD" > to zero, amounts to neglecting magnetisation that starts on the ground state > ensemble and end on the excited state ensemble and vice versa. This will be > a good approximation when pA >> p_B. > > @@ -98,11 +98,87 @@ > """ > > # Python module imports. > -from numpy import arccosh, cos, cosh, sqrt > +import numpy > +from math import cos,sin,atan2 > > > -def r2eff_B14(r20a=None, r20b=None, pA=None, dw=None, kex=None, > cpmg_frqs=None, back_calc=None, num_points=None): > - """ > +def r2eff_B14(r20a=None, r20b=None, pA=None, dw=None, kex=None, power=None, > relax_time=None, tcp=None, back_calc=None, num_points=None): > + """Calculate the R2eff values for the CR72 model. > + > + See the module docstring for details. > + > + > + @keyword r20a: The R20 parameter value of state A (R2 with no > exchange). > + @type r20a: float > + @keyword r20b: The R20 parameter value of state B (R2 with no > exchange). > + @type r20b: float > + @keyword pA: The population of state A. > + @type pA: float > + @keyword dw: The chemical exchange difference between states > A and B in rad/s. > + @type dw: float > + @keyword kex: The kex parameter value (the exchange rate in > rad/s). > + @type kex: float > + @keyword power: The matrix exponential power array. The number > of CPMG blocks. > + @type power: numpy int16, rank-1 array > + @keyword relax_time: The total relaxation time period (in seconds). > + @type relax_time: float > + @keyword tcp: The tau_CPMG times (1 / 4.nu1). > + @type tcp: numpy rank-1 float array > + @keyword back_calc: The array for holding the back calculated R2eff > values. Each element corresponds to one of the CPMG nu1 frequencies. > + @type back_calc: numpy rank-1 float array > + @keyword num_points: The number of points on the dispersion curve, > equal to the length of the cpmg_frqs and back_calc arguments. > + @type num_points: int > """ > > - return None > + # Conversion from relax parameters, to the exact code of Baldwin. > + pb = 1 - pA > + R2e = r20b > + R2g = r20a > + Trel = relax_time > + ncyc = power > + > + ######################################################################### > + ##### Baldwins code. > + ######################################################################### > + pa=(1-pb) > + keg=kex*(1-pb) > + kge=kex*pb > + deltaR2=R2e-R2g > + # This is not used > + #nu_cpmg=ncyc/Trel > + #tcp=Trel/(4.0*ncyc) #time for one free precession element > + > + ######################################################################### > + #get the real and imaginary components of the exchange induced shift > + g1=2*dw*(deltaR2+keg-kge) #same as carver richards zeta > + g2=(deltaR2+keg-kge)**2.0+4*keg*kge-dw**2 #same as carver richards psi > + g3=cos(0.5*atan2(g1,g2))*(g1**2.0+g2**2.0)**(1/4.0) #trig faster than > square roots > + g4=sin(0.5*atan2(g1,g2))*(g1**2.0+g2**2.0)**(1/4.0) #trig faster than > square roots > + ######################################################################### > + #time independent factors > + N=complex(kge+g3-kge,g4) #N=oG+oE > + NNc=(g3**2.+g4**2.) > + f0=(dw**2.+g3**2.)/(NNc) #f0 > + f2=(dw**2.-g4**2.)/(NNc) #f2 > + #t1=(-dw+g4)*(complex(-dw,-g3))/(NNc) #t1 > + t2=(dw+g4)*(complex(dw,-g3))/(NNc) #t2 > + t1pt2=complex(2*dw**2.,-g1)/(NNc) #t1+t2 > + oGt2=complex((deltaR2+keg-kge-g3),(dw-g4))*t2 #-2*oG*t2 > + Rpre=(R2g+R2e+kex)/2.0 #-1/Trel*log(LpreDyn) > + E0= 2.0*tcp*g3 #derived from relaxation #E0=-2.0*tcp*(f00R-f11R) > + E2= 2.0*tcp*g4 #derived from chemical shifts > #E2=complex(0,-2.0*tcp*(f00I-f11I)) > + E1=(complex(g3,-g4))*tcp #mixed term (complex) (E0-iE2)/2 > + ex0b=(f0*numpy.cosh(E0)-f2*numpy.cos(E2)) #real > + ex0c=(f0*numpy.sinh(E0)-f2*numpy.sin(E2)*complex(0,1.)) #complex > + ex1c=(numpy.sinh(E1)) #complex > + v3=numpy.sqrt(ex0b**2.-1) #exact result for v2v3 > + y=numpy.power((ex0b-v3)/(ex0b+v3),ncyc) > + v2pPdN=(( complex(deltaR2+kex,dw) )*ex0c+(-oGt2-kge*t1pt2)*2*ex1c) > #off diagonal common factor. sinh fuctions > + Tog=(((1+y)/2+(1-y)/(2*v3)*(v2pPdN)/N)) > + > Minty=Rpre-ncyc/(Trel)*numpy.arccosh((ex0b).real)-1/Trel*numpy.log((Tog.real)) > #estimate R2eff > + > + # Loop over the time points, back calculating the R2eff values. > + for i in range(num_points): > + > + # The full formula. > + back_calc[i] = Minty[i] > > > _______________________________________________ > relax (http://www.nmr-relax.com) > > This is the relax-commits mailing list > [email protected] > > 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-commits _______________________________________________ relax (http://www.nmr-relax.com) This is the relax-devel mailing list [email protected] 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
