# [gmx-developers] Bug in g_hbond

Michael Patra patra at lorentz.leidenuniv.nl
Thu Mar 3 11:53:22 CET 2005

```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 = 7.91944504
(gdb) print box
\$2 = 7.91944551
(gdb) print ngrid
\$3 = 18
(gdb) print grididx
\$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

```

More information about the gromacs.org_gmx-developers mailing list