colvars icon indicating copy to clipboard operation
colvars copied to clipboard

Minimum image code in NAMD/VMD is not robust to strongly tilted PBCs

Open jhenin opened this issue 6 years ago • 7 comments

The minimum image code found in position_distance() assumes angles close to 90°. In one example, it fails with an angle at 72° and a distance a little smaller that the half box size, but works if the angle is set to 73°.

This could be handled as an error condition, or with more sophisticated PBC code.

jhenin avatar Jun 26 '18 14:06 jhenin

To elaborate (also following private chat), the NAMD function uses the Lattice::delta() method from NAMD, i.e. the error arises from NAMD itself. The function colvarproxy::position_distance() referenced here is only used by VMD right now.

giacomofiorin avatar Jun 26 '18 14:06 giacomofiorin

In-person conversation JH & GF: need to bring this up with NAMD developers.

jhenin avatar Feb 28 '19 15:02 jhenin

Here is a minimalist example. Pathological behavior happens around the half box size. Here, a distance that is exactly half the box length in y becomes 0. PBC_bug_example.zip

jhenin avatar Mar 27 '19 17:03 jhenin

Comment from Jim Phillips:

The Lattice::delta() function is meant for things like bonded terms. If you want the true nearest image you want a variant of Lattice::wrap_nearest_delta().

jhenin avatar Apr 04 '19 16:04 jhenin

That function adds a loop over the 26 nearest unit cells, so it may work for tilt factors a bit over 1/2 (60°), but probably not much higher. @jhenin How much tilt factor were you worried about here?

Note that the additional loop will add a bit of overhead, so a variant will be more useful than the general function (probably a template, with a flag indicating triclinic PBCs as one of its arguments).

giacomofiorin avatar Apr 04 '19 17:04 giacomofiorin

I suppose the question is: how much tilt factor should Colvars support? The case that initially brought the issue to my attention was a truncated octahedral box. I think it is reasonable to support any case that is covered by a search over neighboring unit cells.

jhenin avatar Apr 04 '19 17:04 jhenin

Agreed. There should be an error thrown if the tilt factor exceeds 1 in absolute value.

giacomofiorin avatar Apr 04 '19 18:04 giacomofiorin