[gmx-developers] Bug in g_hbond

[Michael Patra]

Dear all,

I want to report what I believe to be a bug in g_hbond, leading to a
segmentation violation in build_grid(). The relevant parts of the source
code are as follows:

      invdelta[m]=ngrid[m]/box[m][m];

      /* put atom in the box */
      if( x[ad[i]][m] < 0 ) 
        rvec_inc(x[ad[i]],box[m]);
      if( x[ad[i]][m] >= box[m][m] ) 
        rvec_dec(x[ad[i]],box[m]);
      /* determine grid index of atom */
      grididx[m]=x[ad[i]][m]*invdelta[m];
	    
      newgrid=&(grid[grididx[ZZ]][grididx[YY]][grididx[XX]].d[gr]);

The particles are put onto a grid containing ngrid[] elements. The conditional
statements make sure that coordinate x[] of the particle is inside the box. At
first sight, this should give a grid index 0<=grididx[]<ngrid[] such that the
indices into the array grid[] are well-defined.

Unfortunately this is not the case due to the finite precision of invdelta when
single precision floating-point numbers are used. If x[] is close to box[], the
grid index actually can be too large:

(gdb) print x[ad[i]][1]
$1 = 7.91944504
(gdb) print box[1][1]
$2 = 7.91944551
(gdb) print ngrid[1]
$3 = 18
(gdb) print grididx[1]
$4 = 18

The solution to the problem is obvious:

        grididx[m]=x[ad[i]][m]/box[m][m]*ngrid[m];

In addition, I wonder why x[] can be outside the box by at most a single
box size - especially if trjconv -pbc nojump is used. Perhaps there is some
pre-processing of the input data done somewhere in one of the library functions
(I didn't check it). If not, perhaps one should replace the if-statements
by while-statements.

      while( x[ad[i]][m] < 0 ) 
        rvec_inc(x[ad[i]],box[m]);
      while( x[ad[i]][m] >= box[m][m] ) 
        rvec_dec(x[ad[i]],box[m]);

Greetings,

- Michael Patra