[gmx-developers] g_bond can mismeasure distances

[Julius Su]

Hi everyone,

I ran across an issue that can cause the GROMACS g_bond utility to 
mismeasure distances, when the following conditions are satisfied:

(1) The calculation being analyzed was performed without periodic 
boundary conditions.
(2) Average distance calculation is turned on (-averdist) -- this is the 
default.
(3) Atoms are separated by greater than half the distance of the 
bounding box.

In this case, the distances will be measured as if periodic bounds were 
in place, using the minimum image convention. This can lead to 
abnormally low bond lengths, particularly for small molecules.

The issue appears caused by the following piece of code in gmx_bond.c:

   if (bAverDist)
     fdist = opt2fn("-d",NFILE,fnm);
   else {
     fdist = opt2fn_null("-d",NFILE,fnm);
     if (fdist)
     {
       read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&x,NULL,box,
                     FALSE);
     }
   }

which causes periodic bound information to be read in by read_tps_conf 
when bAverDist is true. When bAverDist is false, the read_tps_conf is 
skipped entirely, and ePBC defaults to -1. This causes the periodic 
bounds to be estimated from the stored box size, so that it will be 
inappropriately turned on -- even if previously explicitly turned off.

I propose that the code be changed to:

   if (bAverDist)
     fdist = opt2fn("-d",NFILE,fnm);
   else {
     fdist = opt2fn_null("-d",NFILE,fnm);

   if (fdist)
   {
       read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&x,NULL,box,
                     FALSE);
   }

I am not 100% certain this would not have some detrimental effect 
further downstream. Could one of the developers verify?

Thanks,

Julius