[gmx-developers] g_bond can mismeasure distances

[Mark Abraham]

On 11/11/2011 4:57 AM, Julius Su wrote:
> 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?

The trajectory formats all allow a periodic box to be specified, but 
IIRC none of them permit an indication of whether that box is 
meaningful. So if the trajectory was not generated by a simulation (e.g. 
trjconv -pbc nojump) the box information can be misleading.

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

The possibility of broken molecules after a simulation is a side-effect 
of domain decomposition parallelism, which could have been of a periodic 
or non-periodic system. Or they could arise from trjconv -pbc res/atom.

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

I think the limitations of the multiple file formats and the things 
people want to do make (1) an unachievable goal.

Mark