Re: CST branch -- April 12, 2011 - 15:28

Hi,

it seems possible, but I think that it means that there will
be for each type of the interaction an unique (d,d2)ri_comps function.
It will be defined in the setup_equation (mf.py, class Mf). Currently,
the (d,d2)ri_comps function is selected just based on the fact whether
the values of physical quantities (internuclei distance, csa) are
optimized or not. Now it will be selected again according to the same
criteria and moreover according the type of interaction. This ensure
that inside the function (d,d2)ri_comps the function
create_(d,d2)ri_comps will be called with appropriate parameters. 

Generally, it is probably not very convenient to optimize all
internuclei distances, but only the most crucial one (to the bonded
hydrogen) would be sufficient. However, the selection should be upon
user. The optimalization of CST will be also quite difficult, because 5
parameters would be needed to optimize (3 euler angles + 2 csa values -
corresponding to two pseudo chemical shielding tensors which stand for
asymmetric chemical shielding tensor - hence I always mention csa1 and
csa2 as well as the cross-term, because they are not independent). So
far I would consider only the optimalization of the eigenvalues of the
tensors, the optimalization of the orientation might be added later.

Best regards,
Pavel

On 22 March 2011 12:23, Edward d'Auvergne <edward@xxxxxxxxxxxxx> wrote:

>> > At this point we would like to address a related question. Currently the
>> > calculation of physical constant is done in a several steps. First, the
>> > physical constant is calculated and the value is stored in the
>> > data.dip_const_func or data.csa_const_func (grad, hess). Then, when the
>> > relaxation rates are calculated, the physical constant is recalculated
>> > by
>> > the function create_dip_func or create_csa_func (grad, hess) (method
>> > setup_equations in class Mf, maths_fns/mf.py).
>> >
>> >             comp_dip_const_func(data, data.bond_length)
>> >             comp_csa_const_func(data, data.csa)
>> >             for i in xrange(data.num_ri):
>> >                 data.dip_comps_func[i] = data.dip_const_func
>> >                 if data.create_dip_func[i]:
>> >                     data.dip_comps_func[i] =
>> > data.create_dip_func[i](data.dip_const_func)
>> >                 if data.create_csa_func[i]:
>> >                     data.csa_comps_func[i] =
>> > data.create_csa_func[i](data.csa_const_func[data.remap_table[i]])
>> >
>> > There is one exception, the dipolar physical constant is not
>> > recalculated in
>> > the case of calculation R1 relaxation rate, because the function
>> > create_dip_func does not exist in this case. We do not see a reason for
>> > such
>> > a recalculation.
>>
>> The reason is because of the m10 to m39 models built into relax.  I
>> have made it possible to optimise the bond length and CSA information.
>>  However these models are not stable with the current relaxation data.
>>  I do plan on working with these in the future though, so it would be
>> useful to keep them.  Note that for models m0 to m9, the
>> data.create_dip_func[i] and data.create_csa_func[i] function pointers
>> are set to None.  Therefore for normal model-free analysis the
>> constant is not recalculated.
>>
>>
>> > It seems better to us to just change the coefficient in the
>> > functions comp_r1_dip_jw, comp_r2_dip_jw, comp_r1_csa_jw, comp_r2_csa_jw
>> > (maths_fns/ri_comps.py). I guess, that this design was dedicated to
>> > avoid
>> > multiple calculation of the same interaction constant for each measured
>> > relaxation rate. We would suggest to reach the same effect by this
>> > construction:
>> >
>> >             for i in xrange(data.num_ri):
>> >                 if data.const_func[0]:
>> >                     data.const_func[i] = data.const_func[0]
>> >                 else
>> >                     data.create_const_func(data)
>>
>> For models m10 to m39, this construct will not work.  The constants
>> are already pre-calculated for models m0 to m9 so this is not needed.
>>
>>
>> > Note, the comp_dip_const_func and comp_csa_const_func should be change
>> > so
>> > that, it is possible to call them just with the argument data
>> > (maths_fns/ri_comps.py). Instead of:
>> >
>> >
>> > def comp_dip_const_func(data, bond_length):
>> >     """Calculate the dipolar constant.
>> >
>> >    ...
>> >
>> >     if bond_length == 0.0:
>> >         data.dip_const_func = 1e99
>> >     else:
>> >         data.dip_const_func = 0.25 * data.dip_const_fixed *
>> > bond_length**-6
>> >
>> >
>> > It should look like:
>> >
>> > def comp_dip_const_func(data):
>> >     """Calculate the dipolar constant.
>> >
>> >    ...
>> >
>> >     if data.internuclei_distance == 0.0:
>> >         data.const_func = 1e99
>> >     else:
>> >         data.const_func = 0.25 * data.dip_const_fixed *
>> > data.internuclei_distance**-6
>>
>> The bond_length arg was designed so that the bond length could either
>> come from a fixed value supplied by the user or from the parameter
>> vector when bond lengths or CSA values are optimised.  This behaviour
>> might have to be preserved.
>>
>>
>> > data.dip_const_func were renamed to more general data.const_func and
>> > instead
>> > of bond_length the function directly takes  the internuclei distance for
>> > the
>> > current dipole-dipole interaction. The change of data.dip_const_func to
>> > data.const_func later simplify the code design in the
>> > maths_fns/ri_prime.py
>> > . It will be reduced just to a multiplication of constant and the linear
>> > combination of spectral density functions.
>>
>> For models m10 to m39, I'm not sure if this design would work.  Could
>> we redesign this in another way in which these complex models are
>> still functional?
>>
>
> Then, we would suggest to call the function  comp_dip_const_func,
> comp_csa_const_func ... with full set of possible parameters, i.e.
> comp_dip_const_func(data,internuclei_distance,csa1,csa2,rex)
> comp_csa_const_func(data,internuclei_distance,csa1,csa2,rex)
> ...
> If we would call the function with just physically relevant arguments then
> we would have to use yet another condition to decide the type of the
> interaction in the loop where individual interaction contributions are
> calculated. The physical quantities irrelevent for the given interaction
> (for example csa1,csa2,rex for the dipole-dipole interaction) are None from
> the initialization and will not be used by the function anyway.

It would be easier to have a different function for each physical
constant, as this is very specific code.  For example
comp_dip_const_func is only called from functions which require the
dipolar constant.  For the normal model-free models, this is only
called during setup.  So there is no need to pass in the csa1 and csa2
constants, or Rex for this specific function.  For the CSA constants,
do you calculate 2 separate constants?  For example would you use:

comp_dip_const_func(data[i],internuclei_distance)
comp_csa_const_func(data[j],csa1)
comp_csa_const_func(data[k],csa2)

where data[x] in each case is a different interaction?

>> > Moreover, there is an unanswered question about the NOE and the
>> > additional
>> > dipolar interaction. I am not sure if the suggested design is physically
>> > correct, rather not. During the NOE experiment, the protons are
>> > saturated in
>> > order to reach the steady state. Then a complete set cross relaxation
>> > rates
>> > between all interacting spin pairs should be taken into the account, not
>> > only between the spin of interest and all other interacting nuclei. On
>> > the
>> > other hand this is probably beyond the aim of the program relax. What do
>> > you
>> > think about that?
>>
>> This is getting quite complex as you would need to take the
>> cross-correlated relaxation rates between the different relaxation
>> interactions into account, as well as the motion of all spins if they
>> are not directly bonded.  Is this needed for the current work?  Of
>> course anything is accepted into relax, especially if you would like
>> to probe this area (with a paper in mind), but it has to play nicely
>> with the rest of relax and not be a burden on the relax developers to
>> maintain in the future (as well as not make the current number
>> crunching code slower than it already is).  The code would therefore
>> need to be designed in public.  So if you would like to tackle such a
>> task, I would first recommend finishing off the cst branch, and then
>> make a new branch for this work.
>
> Then, I will design the code as I suggested. So, the sigma_noe will be
> calculated separately for all dipole-dipole interactions with the central
> spin (assuming they are isolated spin pairs). The total sigma_noe will be
> calculated as a sum of all individual sigma_noe contributions.

This sounds the most reasonable for a first approximation.  In reality
the cross rates will be important as well.

Cheers,

Edward