[gmx-developers] add a new autocorrelation function.

[Mark Abraham]

Xavier Periole wrote:
> 
> Dears,
> 
> I am trying to add a new autocorrelation function. The idea if to get
> the residence time of water molecules, but it does not really matter.
> 
> I am modifying gmx3.3.1.
> 
> I made my modifications in gmx_analyze.c to introduce a new option that
> calls do_autocorr the following way:
> 
> if(bacNewoption) {
>   do_autocorr(acfile,"Autocorrelation",n,nset,val,dt,
>                 eacNewOption,bAverCorr);
> }
> 
> I made the same modification for another case where no new function
> was needed and it worked pretty good: reproducing C(t) from another
> source. So this is safe.
> 
> in autocorr.c I also make my changes, defined my variables, etc.
> 
> Main change:
> ...
>       if (MODE(eacNormal)) {
>       corr[k] += c1[j]*c1[j+k];
>       }
>       else if (MODE(eacNewOption) {
>       fprintf(stdout,"in there");
>       compute the new correlation function;
>       }
> ...
> 
> Everything goes well, no complain at the compilation.
> g_analyze calls do_autocorr with the new option but when in autocorr.c
> the code in eacNewOption is not executed. The option is recognized.
> 
> I find out that the options are defined in include/gstat.h in the
> following way:
> 
> #define eacNormal (1<<0)
> #define eacCos    (1<<1)
> #define eacVector (1<<2)
> #define eacRcross (1<<3  | eacVector)
> #define eacP0     (1<<4  | eacVector)
> #define eacP1     (1<<5  | eacVector)
> #define eacP2     (1<<6  | eacVector)
> #define eacP3     (1<<7  | eacVector)
> #define eacP4     (1<<8  | eacVector)
> 
> I do not know how to include my eacNewOption in there! I tried several
> things but the new option is not understood.
> 
> Anyone could give me a hand?

I'm presuming you know what the bitwise-left-shift << and bitwise-or | 
operators do.

You can see in autocorr.c that these get used in a macro test 
"MODE(eacVector)", which is defined in that file and expands to "(mode & 
(eacVector)) == (eacVector)". If eacVector contains particular bits set 
(often only one, as here in the definition from gstat.h), then the test 
evaluates as true if and only if mode has the bits defined in eacVector 
set. This is called "masking" and eacVector is a "bitmask". mode can 
also have other bits set, and from the gstat.h definition, you might 
infer that eacRcross and eacP* are particular instances of eacVector. So 
MODE(eacVector) will be true if mode == eacP4, but not if mode == 
eacCos. This kind of thing is a commonly-used strategy in C to have a 
single flag variable that controls the detail of an operation, and 
GROMACS uses it all over the place. Here's some info which might be 
helpful http://de-coded.blogspot.com/2006/05/basic-bit-masks.html

Since mode is an unsigned long, there are many more bits to play around 
with, and you can try

#define eacNewOption (1<<9)

and see if it works. You may or may not want the "| eacVector", but 
you'll have to judge that from understanding the detail of the code, 
which I don't.

Cheers,

Mark