Re: [sr #3155] An R1rho expression for a spin in chemical exchange between two sites with unequal transverse relaxation ratesRe: [sr #3155] An R1rho expression for a spin in chemical exchange between two sites with unequal transverse relaxation rates|
On 30/04/2014 16:16, Edward d'Auvergne
wrote:
Hi Andy, Do you have code for this R1rho model? It may help to have the code, either to use it directly or for comparison and test suite purposes. In relax, the implementation is in Python. Here some very unoptimised code. Note that you have to adjust the offset (see below). J. Biol. NMR (2013), 55:211-218. In the paper is a formula that predicts accurately when it will work and when it will not (Page 214, 4 lines from the bottom). This formula will break when the relevant eigenvalue is no longer the smallest. As ever, this treatment proves that the R1rho is resilient to differences in relaxation, but to analyse data, use at minimum the 6x6 matrix. But as numpy is used, the different between Python and C is minimal. The fine details of the implementation itself make more of a difference for speed than the choice of language. Debatable. Sure a decent workup in numpy gets you much faster, but you'll notice even the ardent supports of numpy concede that C whoops it's ass. Ie optimised numpy will be slower than C. I say that as someone that does almost everything in python. For the spin-lock offset, it is input and stored as ppm units. The ppm to rad.s^-1 conversion code is visible at http://www.nmr-relax.com/api/3.1/lib.nmr-pysrc.html#frequency_to_rad_per_s. Then the w_eff values are calculated from this in the code at http://www.nmr-relax.com/api/3.1/specific_analyses.relax_disp.disp_data-pysrc.html#return_offset_data. Can you see a bug there? Are you talking about the +/- 0.5*delta_omega factor for slow exchange? Or the Korzhnev 2005 correction for constant time R1rho experiments (in the TODO list at http://www.nmr-relax.com/manual/do_dispersion_features_yet_be_implemented.html, hence it is not implemented yet). There's another trick that doesn't seem to be there. Note that the observed peak position (O_obs) is not the same as the ground state peak position (O_G). Exchange of course moves our peaks around. Experimentally you'll see your peak, and work out where to put the carrier using this as a reference (eg, 100Hz north of peak BLAH). So the apparent measured offset from the data is O_obs not O_G (which I think you assume here). In the analysis, you can account for this simply by subtracting the exchange induced shift from the observed resonance position (taking care to get the sign right). For doing things like determining signs using R1rho eg: J. Biol. NMR (2012) 53(1): 1 You very much need this for example, if you have small deltaOmegas, and/or small offsets and you want to get things right. Again, if you adjust the offset used by the equations by the exchange induced shift, and use this value in the formulas, you'll get an answer that best reflects the data coming off the spectrometer. CATIA is useful for the handling of CPMG off-resonace effects, but it is limited to the models it supports. It cannot handle MMQ data (this is a shortcut for the combination of SQ, ZQ, DQ, and MQ data for both the proton and heteronumcleus). Depends on the version. Mine can do MQ data. You'll notice that in the modern structures from the Kay group almost all use SQ data. This is much easier to handle. Dimitri still likes getting MQ data, but I think he's in the minority in that. CATIA also does not handle R1rho data or models. Anyway, relax has a user function called relax_disp.catia_input which can simplify the life of users by creating the required CATIA input files. There is even a relax_disp.catia_execute user function to allow CATIA to be run from within relax. That's actually a great feature I think. Good stuff! How are you dealing with e.g. selecting the right cpmg sequence and basis sets in the input scripts generated by your code? As detailed in my last email, neglecting this stuff will give you wrong answers. These things have been standard in Flemming's CATIA since its inception. This software was used for the excited state structures you might have seen from the Kay lab recent that have encouraged a great many people to try CPMG. So if people want to emulate that, probably the responsible thing to do is to encourage them to do the analysis rigorously.If a relax_disp.catia_extract user function is added, then CATIA could be used as a replacement optimisation engine in relax and the user will only need to have an executable CATIA binary in their system path. A similar set of user functions exist for Art Palmer's CPMGFit, NESSY and ShereKhan. For relax, the idea is to add the off-resonance support to a few of the numeric models (not all as the Tollinger & Skrynnikov Maple expansion cannot support it). This is in the todo list (http://www.nmr-relax.com/manual/do_dispersion_features_yet_be_implemented.html). As no other dispersion software supports the scalar coupling and spin flips, as far as I am aware, there are no plans or pressure to add this. Best, Andy. Cheers, Edward On 30 April 2014 12:17, Andrew Baldwin <andrew.baldwin@xxxxxxxxxxxxx> wrote:R1rho: Would you prefer python or C? Software seems a bit of a stretch. People code up formulas and use them. Note with R1rho there is a trick that is easily missed. I wondered when you specify an offset to the program in Hz, how precisely your code deals with that? Ie, how precisely do you use that number to get deltaA and OmegaA? Probably you're doing it right, but I didn't see any mention of the trick, which case me to be concerned. Btw, I think your project is a good one! Particularly the intrinsic analysis parts. If you incorporate a few more bits into the CPMG as i've ranted previously, then I think you'll have made something genuinely useful. Catia has everything in, but is a dog to use. Those in the Kay group with the source code get it, but I think no-one outside that world has used it. Which is a shame, as the papers with excited state structures etc that get people interested in CPMG, come from that approach. Best, Andy. On 30/04/2014 11:09, Edward d'Auvergne wrote:Hi Troels, Maybe you should add this BK13 model, as well as B14 (http://wiki.nmr-relax.com/B14) to the TODO section of the dispersion chapter of the relax manual? It could also be added to the table in the docs/latex/dispersion_software.tex file for completeness. I am unaware of any software which implements these models from Andy to add ticks to that table (table 10.3 of http://download.gna.org/relax/manual/relax.pdf). Cheers, Edward On 30 April 2014 00:59, Troels E. Linnet <NO-REPLY.INVALID-ADDRESS@xxxxxxx> wrote:URL: <http://gna.org/support/?3155> Summary: An R1rho _expression_ for a spin in chemical exchange between two sites with unequal transverse relaxation rates Project: relax Submitted by: tlinnet Submitted on: Tue 29 Apr 2014 10:59:47 PM UTC Category: None Priority: 5 - Normal Severity: 3 - Normal Status: Postponed Assigned to: None Originator Email: Open/Closed: Open Discussion Lock: Any Operating System: None _______________________________________________________ Details: In the interests of formula collection, this might be useful for R1rho: An R1rho _expression_ for a spin in chemical exchange between two sites with unequal transverse relaxation rates. J. Biol. NMR (2013), 55:211-218 dx.doi.org/10.1007%2Fs10858-012-9694-6 It's applicable in the same regime as TP02 http://wiki.nmr-relax.com/TP02 Except it will be correct when R2E is not equal to R2G, where other R1rho models may trip up in this regime. _______________________________________________________ Reply to this item at: <http://gna.org/support/?3155> _______________________________________________ Message sent via/by Gna! http://gna.org/ _______________________________________________ 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 |
#!/usr/bin/python
################################################################
# Script to calculate offsets using R1rho experiment
#
# A.Baldwin 12th October 2011
# Copyright (c)
# Please don't circulate without permission
import os,sys,math,string,copy
from math import pi, sin, cos,atan,exp,tan
##########################################################################
def GetTheta(gB1,offset): #return the angle in radians
if(gB1==0):
th=0.0
elif(offset==0 and gB1!=0):
th=90.0/180.0*pi
else:
th=atan(gB1/offset);
return th
def CalcR1rho_palmer(gB1,offset,dw,kex,pb,R1,R2):
th=GetTheta(gB1,offset)
Rex=(pb*dw**2.0*kex)/((offset+dw)**2.0+gB1**2.0+kex**2.0); # according
to Trott/Palmer for pa >> pb
R1rho=R1+(R2+Rex-R1)*(sin(th))**2.0;
return R1rho
def CalcR1rho_New(gB1,offset,dw,kex,Pb,R1,R2,R2E):
#making offset go in opposite direction for spectrometer
deltaA=-offset
deltaB=-offset+dw
deltaBar=-offset+Pb*dw
theta=atan(gB1/deltaBar)
thetaflip=atan(gB1/deltaA)
# for dodgy approximation:
# deltaBar=deltaA
# theta=thetaflip
OmegaA=(gB1**2.+deltaA**2.)**0.5
OmegaB=(gB1**2.+deltaB**2.)**0.5
OmegaE=(gB1**2.+deltaBar**2.)**0.5
Pa=(1-Pb)
R=R2-R1
dr=R2E-R2
#Palmer original terms (not dependent on relaxation difference)
f1p=+ Pa*Pb*dw**2 #sin
f2p=+ kex**2+gB1**2 + deltaB**2.*deltaA**2./deltaBar**2. #cos
dp =+ kex**2 + OmegaA**2*OmegaB**2/OmegaE**2 #Palmer
denominator
#terms affected by relaxation difference
f1= (OmegaA**2 + kex**2. + Pa*kex*dr )*Pb #sin # BIG TERM (2)
f2= dr*Pa+2*kex+gB1**2./kex #sin
f3= 3*Pb*kex+( 2*Pa*kex+gB1**2./kex + dr + dr*Pb**2.*kex**2./OmegaA**2.
)*1/sin(thetaflip)**2. #fits neither sin nor cos
#Coefficients
R2ex=(f1p*kex+dr*f1) / (dp+dr*f3*sin(theta)**2.)
CR1=(f2p+(f1p+dr*(f3-f2))*tan(theta)**2.) / (dp+dr*f3*sin(theta)**2.)
CR2=(dp/sin(theta)**2.-f2p/tan(theta)**2.-f1p+dr*f2) /
(dp+dr*f3*sin(theta)**2.)
#composite expression
R1rhoNew=CR1*R1*math.cos(theta)**2.+(CR2*R2+R2ex)*math.sin(theta)**2.
#R2ex=(dr*Pb) / (1+dr/kex) #limit when w1 is large
return R1rhoNew
theta=atan(gB1/deltaA) #using theta=thetaflip
OmegaA=(gB1**2.+deltaA**2.)**0.5
OmegaB=(gB1**2.+deltaB**2.)**0.5
Pa=(1-Pb)
R=R2-R1
dr=R2E-R2
####approx###########
#setting theta=thetaflip, PG>>PE, kex>>R2-R1
#dp =+ kex**2 +OmegaB**2. #Palmer denominator
#f2= 2*kex+gB1**2./kex #sin
#f3= 3*Pb*kex+( 2*Pa*kex+gB1**2./kex )*1/math.sin(theta)**2. #fits
neither sin nor cos
#
#R2exap=Pa*Pb*dw**2.*kex/dp -
dr*Pb*(Pa*dw**2.*kex**2.*(Pb*(1.-3*cos(theta)**2.)+ 2.+gB1**2./kex**2. ) -
(OmegaA**2.+kex**2.)*(kex**2.+OmegaB**2.))/( (OmegaB**2.+kex**2.)**2.)
#dropping (3cos^2theta-1 )* pb term
R2exap=(Pa*Pb*dw**2.*kex)/(kex**2.+OmegaB**2.)*(1+dr/kex*((kex**2.+OmegaA**2.)/(Pa*dw**2.)-(2*kex**2.+gB1**2.)/(kex**2.+OmegaB**2.)))
CR1ap=1. + dr*Pb*kex*math.tan(theta)**2.*(1.- 3.*cos(theta)**2. ) /
(kex**2.+OmegaB**2.)
CR2ap=1. - dr*Pb*kex*(1-3*cos(theta)**2.) / (kex**2+OmegaB**2.)
#composite expression
R1rhoNew=CR1ap*R1*math.cos(theta)**2.+(CR2ap*R2+R2exap)*math.sin(theta)**2.
#make approximations - gB1=offset and dw^2=offset^2+gB1^2
#print 'OmegaE',((fabs(offset)+dw)**2.+gB1**2.)**0.5,(2+2**0.5)**0.5*dw
#print 'OmegaB',OmegaB,((2-2**0.5)**0.5)*dw
#print 'OmegaA',(offset**2.+gB1**2.)**0.5,dw
#print 'Offset',offset,dw/(2**0.5)
#print 'gB1',gB1,dw/(2**0.5)
#OmegaB=(2-2**0.5)**0.5*dw #if we're going with the negative definition
of offset for omegaB
#OmegaA=dw
#gB1=dw/2**0.5
return R1rhoNew