[gmx-users] Small bug with -pi in g_mindist?

[Reid Van Lehn]

Hi users,

I was experimenting with using a hexagonal unit cell for a lipid membrane
system by building a box with vectors A A C (i.e. two equivalent vectors in
the xy plane and distinct z vector) and angles 90 90 60, which is correctly
represented as hexagonal after converting with trjconv. I equilibrated this
box and a system with the same number of lipids and water molecules but in
a rectangular box. I confirmed that after equilibration the area per lipid
for each system is the same, and after  visualization confirmed that the
hexagonal system and rectangular system occupied about the same area (of
course subject to fluctuations).

For a hexagon and square of equivalent area, the minimum distance between
periodic images should be 1/sqrt( sin 60 degrees) = ~1.075 times larger for
the hexagon than in the square case if I worked out the geometry correctly.
To test to make sure I had the correct new minimum distance, I ran
g_mindist -pi with both a single atom and a single water molecule from the
simulation box after resizing the Z axis to a large value, restricting the
minimum distance between periodic images to only the xy plane.
Surprisingly, the result came out as smaller than the equivalent minimum
distance in the rectangular box, and was equal to the box vector B. Since
in GMX box vectors are stored in the .gro file in a 3x3 matrix, the
box[YY][YY] vector for an ab angle of 60 degrees was (correctly) equal to
sin 60 * the A vector. However, the correct periodic distance should have
been the A  vector, which again was correctly ~1.075 * the box vector in
the rectangular box by comparing the .gro files.

I believe that this is a small bug in g_mindist, and found the source: on
lines 71 and 92 of gmx_mindist.c (version 4.6), the initial minimum
distance is set based on the minimum of the box[XX][XX], box[YY][YY], and
box[ZZ][ZZ] vectors. I think this is fine for rectangular boxes, but fails
for triclinic boxes with box angles differing from 90 degrees. In my
particular case of a hexagonal xy plane, the box[YY][YY] vector is shorter
than the box[XX][XX] vector by a factor of sin 60, but this is not actual
the shortest distance. Explicitly printing out the distances calculated
between periodic images in the current version of the code for both a
single atom and a water molecule confirms that the minimum distance is the
box[XX][XX] vector.

A fix for this was substituting norm(box[0]), norm(box[1]), and
norm(box[2]) for the box[XX][XX] etc. vectors in line 71, as then the
minimum was properly set.

I realize this is a bug that is unlikely to affect many users, but given
the prevalence of non-rectangular boxes and the observation of previous
complaints about mindist in the user list, I thought it would be good to
report.

Please let me know if I made an error anywhere!

Best,
Reid