[gmx-developers] g_bond can mismeasure distances

[Mark Abraham]

On 9/11/2011 4:01 PM, Julius Su wrote:
> 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 think your restatement of the logic of this code with respect to 
bAverDist isn't right, but there is certainly a relevant issue to be 
fixed here. See http://redmine.gromacs.org/issues/834.

>
> 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?

I think that is sound, but now g_bond will require -s, and that file 
should be chosen to have a PBC treatment consistent with how your 
analysis is intended to function on that trajectory. There are some 
details there I've never tried to understand, however.

Mark