I tried to run the test with different models. relax -s Relax_disp.test_baldwin_synthetic --time -d
set_grid_r20_from_min_r2eff=True GRID=21 CR72 full http://www.nmr-relax.com/api/3.1/lib.dispersion.cr72-module.html r2a (SQ CPMG - 200.00000000 MHz) 2.856395585456784 r2b (SQ CPMG - 200.00000000 MHz) 4.906962606245907 dw 1.973207453768875 pA 0.991890196690937 kex 1144.199191689991267 time 3.698s NS CPMG 2-site star full http://www.nmr-relax.com/api/3.1/lib.dispersion.ns_cpmg_2site_star-module.html http://wiki.nmr-relax.com/NS_CPMG_2-site_star_full Optimised parameter values: r2a (SQ CPMG - 200.00000000 MHz) 2.645120395040764 r2b (SQ CPMG - 200.00000000 MHz) 2.942150337866721 dw 0.488796499895965 pA 0.912119729304970 kex 3412.892868042056307 time 50.86 s NS CPMG 2-site 3D full http://www.nmr-relax.com/api/3.1/lib.dispersion.ns_cpmg_2site_3d-module.html http://wiki.nmr-relax.com/NS_CPMG_2-site_3D_full Optimised parameter values: r2a (SQ CPMG - 200.00000000 MHz) 1.999959684693684 r2b (SQ CPMG - 200.00000000 MHz) 100.004099440950242 dw 1.999999932481509 pA 0.989999918707243 kex 999.995874168210662 time 134.426s B14 http://wiki.nmr-relax.com/B14 Optimised parameter values: r2a (SQ CPMG - 200.00000000 MHz) 1.999958942472283 r2b (SQ CPMG - 200.00000000 MHz) 100.004175006368143 dw 1.999999933684466 pA 0.989999917242163 kex 999.995802494615418 time 6.001s NS CPMG 2-site expanded http://www.nmr-relax.com/api/3.1/lib.dispersion.ns_cpmg_2site_expanded-module.html http://wiki.nmr-relax.com/NS_CPMG_2-site_expanded r2 (SQ CPMG - 200.00000000 MHz) 2.873100105237559 dw 1.998666628587453 pA 0.991663011056017 kex 1094.588041531584395 time 8.735s 2014-05-03 22:27 GMT+02:00 Troels Emtekær Linnet <[email protected]>: > Dear Andrew. > > I am happy to announce, that I have now implemented your model in relax. > > At this step, I created sample data R2eff data with your provided > script, from the spin dynamic parameters from your paper: > Dynamic parameters > 1) kex=1000. rad/s > 2) pb=0.01 > 3) dw_ppm=2.0 > 4) R2g= r2a = 2.0 rad/s > 5) R2e=r2b = 100. rad/s. > > Experiment settings > 1) isotope('1H', spin_id='@H') > 2) ncycs = [2, 4, 8, 10, 20, 40, 500] > 3) sfrq= 200. * 1E6 # 1H Larmor frequency. > 4) Trelax=0.04 # Total time of CPMG block. > > When loading the generated R2eff data, relax will now do grid search, > minimise, and find the > dynamic parameters. > > This is implemented with the systemtest: > relax -s Relax_disp.test_baldwin_synthetic -d > > (Look for def test_baldwin_synthetic(self): ) > http://svn.gna.org/viewcvs/*checkout*/relax/trunk/test_suite/system_tests/relax_disp.py?revision=HEAD > > An example graph for the test data is here: > http://wiki.nmr-relax.com/B14#Example_graph_of_test_data > > The main code for the implementation is: > > The library code. > > http://svn.gna.org/viewcvs/*checkout*/relax/trunk/lib/dispersion/b14.py?revision=HEAD > > The target function. The func_B14(). > > http://svn.gna.org/viewcvs/*checkout*/relax/trunk/target_functions/relax_disp.py?revision=HEAD > > All commits that created the code, approx 20-30 commits, with > corrections, can be found on the mailing lists with this search. > http://search.gmane.org/?query=https%3A%2F%2Fgna.org%2Fsupport%2F%3F3154&author=&group=gmane.science.nmr.relax.scm&sort=revdate&DEFAULTOP=and&TOPDOC=10&xP=https%09gna%09org%09support%093154&xFILTERS=Gscience.nmr.relax.scm---A > > This follows the tutorial at: > http://wiki.nmr-relax.com/Tutorial_for_adding_relaxation_dispersion_models_to_relax > > Step 1-9 is completed. > > Now follows step 10, which is debugging. > 1) Optimisation of code > 2) Prevent trigonometric/math functions to be overloaded. Currently > the Grid search makes these complaints, for some edge parameter > values. > > RelaxWarning: invalid value encountered in divide > RelaxWarning: invalid value encountered in sqrt > RelaxWarning: invalid value encountered in arccosh > RelaxWarning: invalid value encountered in log > > A good test case would also be the creation of R2eff data with the > already implemented numerical models, and the fit with B14. > > Best > Troels > > > > > > > > > 2014-05-03 14:54 GMT+02:00 Troels Emtekær Linnet <[email protected]>: >> Dear Andrew Baldwin. >> >> I am in the process of implementing your code in relax. >> >> I am getting myself into a terrible mess by trying to compare >> equations in papers, code and conversions. >> >> But I hope you maybe have a little time to clear something up? >> >> I am looking at the current version of your paper: >> dx.doi.org/10.1016/j.jmr.2014.02.023 >> which still is still "early version of the manuscript. The manuscript >> will undergo copyediting, typesetting, and review of the resulting >> proof before it is published in its final form." >> >> 1) Equations are different ? >> In Appendix 1 – recipe for exact calculation of R2,eff >> h1 = 2 dw (dR2 + kEG - kGE) >> >> But when I compare to: >> >> Supplementary Section 4. Relation to Carver Richards equation. Equation 70. >> h1 = zeta = 2 dw alpha_ = 2 dw (dR2 + kGE - kEG) >> >> Is there a typo here? The GE and EG has been swapped. >> >> 2) Missing zeta in equation 70 >> There is missing \zeta instead of z in equation 70, which is probably >> a typesetting error. >> >> 3) Definition of " Delta R2" is opposite convention? >> Near equation 10, 11 you define: delta R2 = R2e - R2g >> if we let e=B and g=A, then delta R2 = R2B - R2A >> >> And in equation 70: >> alpha_ = delta R2 + kGE - kEG >> >> but i CR original work, A8, >> alpha_ = R2A - R2B + kA - kB = - delta R2 + kGE - kEG >> >> So, that doesn't match? >> >> Sorry if these questions are trivial. >> >> But I hope you can help clear them up for me. :-) >> >> Best >> Troels >> >> >> >> >> >> >> 2014-05-01 10:07 GMT+02:00 Andrew Baldwin <[email protected]>: >>> >>> The Carver and Richards code in relax is fast enough, though Troels >>> might have an interest in squeezing out a little more speed. Though >>> it would require different value checking to avoid NaNs, divide by >>> zeros, trig function overflows (which don't occur for math.sqrt), etc. >>> Any changes would have to be heavily documented in the manual, wiki, >>> etc. >>> >>> >>> The prove that the two definitions are equivalent is relatively >>> straightforward. The trig-to-exp command in the brilliant (and free) >>> symbolic program sage might prove helpful in doing that. >>> >>> >>> For the tau_c, tau_cp, tau_cpmg, etc. definitions, comparing the relax >>> code to your script will ensure that the correct definition is used. >>> That's how I've made sure that all other dispersion models are >>> correctly handled - simple comparison to other software. I'm only >>> aware of two definitions though: >>> >>> tau_cpmg = 1 / (4 nu_cpmg) >>> tau_cpmg = 1 / (2 nu_cpmg) >>> >>> >>> tau_cpmg = 1/(nu_cpmg). >>> >>> Carver and Richard's switch definitions half way through the paper >>> somewhere. >>> >>> >>> What is the third? Internally in relax, the 1 / (4 nu_cpmg) >>> definition is used. But the user inputs nu_cpmg. This avoids the >>> user making this mistake as nu_cpmg only has one definition. >>> >>> >>> I've seen people use nu_cpmg defined as the pulse frequency. It's just an >>> error that students make when things aren't clear. I've often seen brave >>> student from a lab that has never tried CPMG before do this. Without someone >>> to tell them that this is wrong, it's not obvious to them that they've made >>> a mistake. I agree with you that this is solved with good documentation. >>> >>> >>> >>> You guys are free to use my code (I don't mind the gnu license) or of course >>> implement from scratch as needed. >>> >>> Cheers! For a valid copyright licensing agreement, you'll need text >>> something along the lines of: >>> >>> "I agree to licence my contributions to the code in the file >>> http://gna.org/support/download.php?file_id=20615 attached to >>> http://gna.org/support/?3154 under the GNU General Public Licence, >>> version three or higher." >>> >>> Feel free to copy and paste. >>> >>> >>> No problem: >>> >>> "I agree to licence my contributions to the code in the file >>> http://gna.org/support/download.php?file_id=20615 attached to >>> http://gna.org/support/?3154 under the GNU General Public Licence, >>> version three or higher." >>> >>> >>> I'd like to note again though that anyone using this formula to fit data, >>> though exact in the case of 2site exchange/inphase magnetisation, evaluated >>> at double floating point precision should not be doing so! Neglecting scalar >>> coupling/off resonance/spin flips/the details of the specific pulse sequence >>> used will lead to avoidable foobars. I do see value in this, as described in >>> the paper, as being a feeder for initial conditions for a more fancy >>> implemenation. But beyond that, 'tis a bad idea. Ideally this should appear >>> in big red letters in your (very nice!) gui when used. >>> >>> In relax, we allow the user to do anything though, via the >>> auto-analysis (hence the GUI), we direct the user along the best path. >>> The default is to use the CR72 and 'NS CPMG 2-site expanded' models >>> (Carver & Richards, and Martin Tollinger and Nikolai Skrynnikov's >>> Maple expansion numeric model). We use the CR72 model result as the >>> starting point for optimisation of the numeric models, allowing a huge >>> speed up in the analysis. The user can also choose to not use the >>> CR72 model results for the final model selection - for determining >>> dispersion curve significance. >>> >>> >>> Here's the supp from my CPMG formula paper (attached). Look at the last >>> figure. Maybe relevant. Nikolai's algorithm blows up sometimes when you >>> evaluate to double float point precision (as you will when you have a >>> version in python or matlab). The advantage of Nicolai's formula, or mine is >>> that they won't fail when Pb starts to creep above a per cent or so. >>> >>> Using the simple formula as a seed for the more complex on is a good idea. >>> The most recent versions of CATIA have something analogous. >>> >>> >>> As for the scalar coupling and spin flips, I am unaware of any >>> dispersion software that handles this. >>> >>> >>> CATIA. Also cpmg_fig I believe. In fact, I think we may have had this >>> discussion before? >>> >>> https://plus.google.com/s/cpmg%20glove >>> >>> If you consider experimental validation a reasonable justification for >>> inclusion of the effects then you might find this interesting: >>> >>> Spin flips are also quite relevant to NH/N (and in fact any spin system). >>> The supplementary section of Flemming and Pramodh go into it here for NH/N >>> http://www.pnas.org/content/105/33/11766.full.pdf >>> >>> And this: >>> JACS (2010) 132: 10992-5 >>> Figure 2: >>> r2 0.97, rmsd 8.0 ppb (no spin flips) >>> r2 0.99, rmsd 5.7 ppb (spin flips). >>> >>> The improvements are larger than experimental uncertainties. >>> >>> When we design these experiments and test them, we need to think about the >>> details. This is in part because Lewis beats us if we don't. You can imagine >>> that it comes as a surprise then when we see people neglecting this. In >>> short, the parameters you extract from fitting data will suffer if the >>> details are not there. In the case of spin flips, the bigger the protein, >>> the bigger the effect. In your code, you have the opportunity to do things >>> properly. This leaves the details behind the scenes, so the naive user >>> doesn't have to think about them. >>> >>> >>> And only Flemming's CATIA >>> handles the CPMG off-resonance effects. This is all explained in the >>> relax manual. I have asked the authors of most dispersion software >>> about this too, just to be sure. I don't know how much of an effect >>> these have though. But one day they may be implemented in relax as >>> well, and then the user can perform the comparison themselves and see >>> if all these claims hold. >>> >>> >>> Myint, W. & Ishima, R. Chemical exchange effects during refocusing pulses in >>> constant-time CPMG relaxation dispersion experiments. J Biomol Nmr 45, >>> 207-216, (2009). >>> >>> also: >>> >>> Bain, A. D., Kumar Anand, C. & Nie, Z. Exact solution of the CPMG pulse >>> sequence with phase variation down the echo train: application to R(2) >>> measurements. J Magn Reson 209, 183- 194, (2011). >>> >>> >>> Or just simulate the off-resonance effects yourself to see what happens. For >>> NH you see the effects clearly for glycines and side chains, if the carrier >>> is in the centre around 120ppm. The problem generally gets worse the higher >>> field you go though this of course depends when you bought your probe. You >>> seem to imply that these effects are almost mythical. I assure you that they >>> come out happily from the Bloch equations. >>> >>> Out of curiosity, from a philosophical perspective, I wonder if you'd agree >>> with this statement: >>> >>> "the expected deviations due to approximations in a model should be lower >>> than the experimental errors on each datapoint." >>> >>> >>> >>> Anyway, the warnings about analytic versers numeric are described in >>> the manual. But your new model which adds a new factor to the CR72 >>> model, just as Dimitry Korzhnev's cpmg_fit software does for his >>> multiple quantum extension of the equation (from his 2004 and 2005 >>> papers), sounds like it removes the major differences between the >>> analytic and numeric results anyway. In any case, I have never seen >>> an analytic result which is more than 10% different in parameter >>> values (kex specifically) compared to the numeric results. I am >>> constantly on the lookout for a real or synthetic data set to add to >>> relax to prove this wrong though. >>> >>> >>> I think there's a misunderstanding in what I mean by numerical modelling. >>> For the 2x2 matrix (basis: I+, ground and excited) from the Bloch McConnell >>> equation, you can manipulate this to get an exact solution. Nikolai's >>> approach does this, though his algorithm can trip up sometimes for >>> relatively exotic parameters when you use doubles (see attached). My formula >>> also does this. I agree with you entirely that in this instance, numerically >>> solving the equations via many matrix exponentials is irrelevant as you'll >>> get an identical result to using a formula. >>> >>> My central thesis here though is that to get an accurate picture of the >>> experiment you need more physics. This means a bigger basis. To have off >>> resonance effects, you need a minimum 6x6 (Ix,Iy,Iz, ground and excited). To >>> include scalar coupling and off resonance, you need a 12x12 >>> (Ix,Iy,Iz,IxHz,IyHz,IzHz, ground and excited). Including R1 means another >>> element and so on. The methyl group, for example, means you need 24. >>> >>> So when we use the term numerical treatment, we generally mean a calculation >>> in a larger basis, as is necessary to take into account the appropriate spin >>> physics. There aren't formulas to deal with these situations. In the 6x6 >>> case for example, you need 6 Eigenvalues, which is going to make life very >>> rough for someone brave enough to attempt a close form solution. The Palmer >>> and Trott trick used in 2002 for R1rho is a smart way of ducking the problem >>> of having 6 Eigenvalues, but for CPMG unfortunately you need several >>> Eigenvalues, not just the smallest. >>> >>> The 2x2 matrix that Nikolai and Martin, Carver Richard's and I analyse does >>> not include scalar coupling, as magnetisation is held in-phase (in addition >>> to neglecting all the other stuff detailed above). So it is a good >>> representation for describing the continuous wave in-phase experiments >>> introduced here (neglecting relaxation effects and off resonance): >>> >>> Vallurupalli, P.; Hansen, D. F.; Stollar, E. J.; Meirovitch, E.; Kay, L. E. >>> Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 18473–18477. >>> >>> and here: >>> >>> Baldwin, A. J.; Hansen, D. F.; Vallurupalli, P.; Kay, L. E. J. Am. Chem. >>> Soc. 2009, 131, 11939–48. >>> >>> But I think are the only two where this formula is directly applicable. Only >>> if you have explicit decoupling during the CPMG period do you satisfy this >>> condition. So this is not the case for all other experiments and certainly >>> not true for those used most commonly. >>> >>> Anyhow. Best of luck with the software. I would recommend that you consider >>> implementing these effects and have a look at some of the references. The >>> physics are fairly complex, but the implementations are relatively >>> straightforward and amount to taking many matrix exponentials. If you do >>> this, I think you'd end up with a solution that really is world-leading. >>> >>> As it stands though, in your position, I would worry that on occasion, users >>> will end up getting slightly wrong parameters out from your code by >>> neglecting these effects. If a user trusts this code then, in turn, they >>> might lead themselves to dodgy biological conclusions. For the time being, >>> I'll stick to forcing my students to code things up themselves. >>> >>> All best wishes, >>> >>> Andy. >>> >>> >>> _______________________________________________ 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
