Alex, welcome to the relax development mailing list. The expansion of
relax to handle both the rhombicity of the chemical shift (CS) tensor
and the non-colinearity of the XH bond and principle axis of the CS
tensor would be much welcome. As demonstrated by your ICMRBS poster,
this should be a significant improvement for those studying RNA and
DNA dynamics. It will probably have applications within protein
systems as well.
The idea behind relax is to incorporate all NMR relaxation data
analysis concepts into the program. I'm also trying to set up the
program to be community driven as I may move on to different areas as
I've now submitted in my PhD thesis. I've set up a comprehensive open
source framework for relax so that the program can never die even if I
give up on it (https://gna.org/projects/relax/).
This will be a long post so I will break it into a few sections.
These include the basics for developing relax, the CS tensor input,
locating homonuclei, optimisation of the new equations, and other
stuff.
1. The basics for developing relax.
Coding in relax shouldn't be that difficult. Essentially Python is
similar in flow to C but there are no declarations, no variable
initialisation, and most of the brackets are dropped. Of all the
programing languages I've played with (that includes Perl, Python, C,
debugging and modifying FORTRAN, Basic, a bit of TCL, a little Scheme,
and a few others), Python would be, by far, the simplest language to
code large programs in. It may not be the most efficient but it is
very scientist friendly (being so basic). If you know how to code in
a proper programming language then mimicking the code already present
in relax to produce your own features should be straight-forward.
The most important coding convention in relax is comments. Just under
half of the text in relax is comments. The purpose of this is to
allow new developers to very quickly understand and learn the relevant
portions of the code. I know this is not normal for computer code but
I prefer self documented code. This self documentation should be of
significant help when making additions to the program. As for
subversion, it should take a few hours to learn the basic commands.
To start with only two commands are needed. The first is to check out
the current sources:
$ svn co svn://svn.gna.org/svn/relax/1.3 relax
The 1.3 line is the best place for these additions. The second is to
update your checked out copy to get all the changes which have been
occurring in the repository:
$ svn up
Another two commands which are useful are
$ svn st
which shows you the files you have modified and
$ svn diff
which shows you the actual changes in diff format. If you would like
to make the changes, then the command
$ svn diff > patch
will create a patch file that can be sent to relax-devel. There are
currently four relax developers (probably more in the future, possibly
including yourself) who can help you along. There is Chris MacRaild,
Andrew Perry, Gary Thompson, and myself. The patch can be checked and
suggestions provided. Once the patch fits well into relax then I or
one of the other developers can add it to the repository. After
learning a few of the basic relax conventions then commit access to
the repository can be granted and you can submit the changes into the
repository yourself. The conventions are very basic and are simply
formatting issues. As for ruining code, the good thing about version
control systems such as subversion is that everything is reversible.
Even if you accidentally delete everything out of the repository, it
can all be recovered in a single command.
2. The CS tensor.
To start with, we should break this into two very distinct components
- minimal user input and then relax handling the rest. Currently the
anisotropic component of the CS tensor, the CSA, is input using a
command such as
relax> value.set('m1', -170 * 1e-6, 'CSA')
The 'value.set()' command is used to set various residue specific
values and the entire CS tensor could be input using this approach.
Alternatively the diffusion tensor initialisation command
'diffusion_tensor.init()' may be a framework which could be used. For
the most complex diffusion tensor, the ellipsoid, 6 parameters are
input - the isotropic component of the tensor (Diso); the anisotropic
component of the tensor (Da); the rhombic component of the tensor
(Dr), and the three Euler angles, alpha, beta, and gamma.
As the anisotropic component of the CS tensor is set using
'value.set()', the rhombic component of the tensor could be input in
the same way. It's value could be specified as 'CSR' and stored in
the residue specific variable 'csr'. This system is also more logical
than specifying the eigenvalues of the CS tensor (the x, y, and z
chemical shift components of the tensor). I do have one question
though. Do we need to input the isotropic chemical shift for your
equations? This would be necessary if the equations use the
eigenvalues rather than the anisotropy and rhombicity. Actually it
would be useful if the equations were represented in text format to a
post to relax-devel. This text could be included in the comments to
the new code.
Another important question is whether the user inputs the CS tensor
orientation (Euler angles, quaternions, etc). Is it necessary that
the user inputs these values or are they dependent on the PDB
coordinates? If we can remove this requirement away from the user and
handle it in the code, it would be better. The less the user does,
the better - this is an important philosophy behind relax. Otherwise
we could again use the 'value.set()' function to set parameters such
as 'CS_alpha', 'CS_beta', and 'CS_gamma'.
The most involved solution would be to create a new user function such
as 'csa.init()' which mimics the 'diffusion_tensor.init()' function.
3. Locating homonuclei.
This could be either very easy or quite complex. If only the
neighbouring bonded heteronuclei need to be taken into account then
these could be easily extracted from the PDB. relax uses the function
'pdb()' to load the PDB file which then uses Scientific python to read
the file. This results in a data structure containing all the PDB
info from which the atoms can be individually selected. Otherwise if
atoms from other residues need to be taken into account, an algorithm
for identifying atoms within a certain radius may need to be coded.
The final question is, what is then done with these spins?
4. The optimisation of the new equations.
This code belongs in the directory 'maths_fns'. The complexity of
these additions would be proportional to the complexity of the
equations. To see how the chi-squared value is calculated in relax
from the model-free parameters and CSA, bond length, and Rex values,
have a look at the start of the Values, gradients, and Hessians
chapter of the relax manual. A clear statement of the problem
including which values are fixed, which are optimised, the equations
used, etc, would be of help in determining the scope of the changes
required.
5. Misc.
You said that the model-free parameters are the same for the CS and
dipolar interactions, does this mean that you assume the same spectral
density J(w) for both relaxation mechanisms? This would simplify
implementation. Also relax handles all types of Brownian rotational
diffusion including diffusion as a sphere, spheroid, and ellipsoid.
The original reference has equations for an ellipsoid and you have
derived the equations for a spheroid. Are there equations for the
sphere (isotropic diffusion)?
The best place to dive in to would be setting the rhombic CS tensor
component in 'value.set()'. This will involve modifying code in the
file 'specific_fns/model_free.py'. If you do get lost in the large
code base, don't hesitate to ask questions in this forum.
The changes required for implementing all the concepts could be quite
disruptive. Hence the changes should occur in the 1.3 line in the
repository (http://svn.gna.org/viewcvs/relax/). The 1.2 line is
stable and is now reserved for minor or non-disruptive changes.
I hope that that answers a few of your questions. Sincerely,
Edward
Re: Improving expressions for the CSA interaction