Tutorial for adding relaxation dispersion models to relax

Tutorial

See mail archive

The following is a tutorial for adding new relaxation dispersion models for either CPMG-type or R1rho-type experiments to relax. This includes both the models based on the analytical, closed-form expressions as well as the models involving numerical solutions of the Bloch-McConnell equations.

The tutorial will follow the example of the addition of the 'M61' model already present within relax, pointing to the relevant commits for reference. To see the commit message and the code changes in colour, click on the links found within these commit messages. This specific case is the Meiboom 1961 analytic model for 2-site fast exchange equation for R1rho-type experiments.

Adding the model to the list

Firstly the model should be added to the lists of the specific_analyses.relax_disp.variables module. The model name is stored in a special variable which will be used throughout relax.

Reference commits

• The M61 model
• The M61 skew model

The relax_disp.select_model user function front end

Reference commit: http://article.gmane.org/gmane.science.nmr.relax.scm/17612

The next step is to add the model, its description, the equations for the analytic models, and all references to the relax_disp.select_model user function front end. When the relaxation dispersion chapter of the relax manual is created (this will be the docs/latex/relax_disp.tex file), then the same description should be added there as well.

The relax library

Reference commit: http://article.gmane.org/gmane.science.nmr.relax.scm/17615

Now the dispersion function needs to be added to the relax library (in the lib.relax_disp package). This should be designed as a simple Python function which takes the dispersion parameters and experimental variables, and calculates the R2eff/R1rho values. The module can contain auxiliary functions for the calculation. Some auxiliary functions, if not specific to relaxation dispersion, may be better placed in other locations within the relax library.

The relaxation dispersion functions in the library currently take as an argument a data structure for the back-calculated R2eff/R1rho values and populate this structure. This design is not essential if the target function, described in the next point, handles the library function appropriately. Just look at the files in lib/dispersion to get an idea of the design used.

The dispersion code in the relax library must be robust. This involves identifying parameter values or combinations which would cause failures in the mathematical operations (numerical issues not present in the mathematics must be considered). Note that parameter values of 0 are common within a grid search. It should be decided if the R2eff/R1rho value should be set to zero, to another value, or to something large (e.g. 1e100). For example:

Divisions - always catch zeros in the denominator with if statements, even if you believe that this will never be encountered. Square roots - make sure that the value inside is always > 0. Trigonometric functions - these should be tested for where they are not defined or where the software implementation can no longer handle certain values. For example try cosh(1000) in Python.

In the reference example, the M61 model code was copied from the LM63 module and modified appropriately.

The target function

Reference commits:

http://article.gmane.org/gmane.science.nmr.relax.scm/17616 http://article.gmane.org/gmane.science.nmr.relax.scm/17660 http://article.gmane.org/gmane.science.nmr.relax.scm/17661

The target function is used in optimisation and is a class method which takes as a single argument the parameter vector. This list is changed by the minimisation algorithm during optimisation. The target function should then return a single floating point number - the chi-squared value.

Again in this example, the code for the M61 is copied from the LM63 model and then modified.

Adding support for the parameters

Reference commit: http://article.gmane.org/gmane.science.nmr.relax.scm/17573

This is needed to enable the model. This example is for the CR72 model implementation as the parameters required for the M61 model match those of the preexisting LM63 model.

The relax_disp.select_model back end

Reference commit: http://article.gmane.org/gmane.science.nmr.relax.scm/17622

Now the back end of the relax_disp.select_model user function for the model can be added. This involved identifying the model and constructing the parameter list.

The test suite

Reference commits:

http://article.gmane.org/gmane.science.nmr.relax.scm/17647 http://article.gmane.org/gmane.science.nmr.relax.scm/17648 http://article.gmane.org/gmane.science.nmr.relax.scm/17649 http://article.gmane.org/gmane.science.nmr.relax.scm/17662 http://article.gmane.org/gmane.science.nmr.relax.scm/17663

This step is normally performed as step number 1. This is the most important part that makes sure that the code not only works now, but will continue working for the entire lifetime of the relax project.

The idea is that synthetic data (here for example as Sparky peak lists) is created for the model and added to the test suite directory test_suite/shared_data/dispersion/. This is then used in a system test to check that the code in relax can reproduce the data. It is very important that the code added to the relax library is not used to create the synthetic data!

Comparing to other software

It can happen that a bug present in the lib.dispersion package code is also replicated in the synthetic data. This is not uncommon. Therefore it is very useful to use other software with the test data from step 7 to see if the original parameters can be found. A good example can be seen in the test_suite/shared_data/dispersion/Hansen which contains Dr. Flemming Hansen's CPMG data (see the README file) and the results from different programs including NESSY, relax, CPMGFit, and ShereKhan. The comparison is in the file 'software_comparison'.

Once the relax code is able to find identical or better results than the dispersion softwares, then the values found in the test suite optimisation can be locked in. The assertEqual() and assertAlmostEqual() methods can be used to only allow the test to pass when the correct values are found.

Debugging

This step should not require an explanation. It goes hand-in-hand with steps 7) and 8).

The auto-analysis

The model variable in specific_analyses.relax_disp.variables needs to be imported into the auto_analyses.relax_disp module. This is then used in the write_results() method to output text files and Grace plots of the parameters. Be sure that the model variable is added to each part of this method corresponding to the parameters of the model.

The GUI

The model needs to also be added to the graphical user interface (GUI). This is in the gui.analyses.auto_relax_disp module. The model variable should first be imported. In the __init__() method, it should be decided if the model should be selected by default or if the user should manually select the model during the analysis. If the former, then it should be added to the ds.relax_gui.analyses[data_index].disp_models list.

For the model to be accessible via the GUI, it must be added to the Disp_model_list_cpmg or Disp_model_list_r1rho model list classes (at the bottom of the module). The model variable should be added to the models list, and the list of parameters added to the params list.