[gmx-developers] g_bond can mismeasure distances

[Julius Su]

On 11/10/2011 5:19 AM, Mark Abraham wrote:
> 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.
Thanks Mark for the clear write-up of the issue. If I understand 
correctly, the PBC can be present in two ways: (1) in the trajectory 
file (.trr) and (2) in the structure file (.tpr), and the two will not 
always agree with each other.

A further question though: under what circumstances would the PBC 
information in the trajectory be incorrect? does this have to do with 
the treatment of broken molecules? Perhaps if the molecules get printed 
as broken, there should be PBC info, but if not, it should revert to 
nonperiodic? I ask as someone unfamiliar with the historical development 
of the code, and the design decisions made along the way.

It seems to me that if the PBC information in (1) could be made to be 
always correct, this would be useful for making sure all the downstream 
analyses were then consistent with the computed trajectory .. (with the 
option of overriding this by specifying a structure file with -s, 
perhaps with a warning that the PBC was being overridden).

Julius

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