[gmx-users] g_density -center issues (GMX version 4.6)

[Reid Van Lehn]

Hi users/developers,

I identified 2 issues in the center_coords function of g_density that make
the flag -center incorrectly center the system. These are in addition to
the issues previously raised by Chris Neale (reported as issue 1168 on
redmine: http://redmine.gromacs.org/issues/1168). I have added these two
issues to the same redmine bug report as they are related to the same
analysis tool/similar issues; following from Chris' previous advice, I also
suggest that users not use the -center flag and instead center their
trajectories prior to running g_density using trjconv. The issues described
below are present in g_density for GMX 4.6 and presumably for earlier/later
versions as well; I apologize if these have been updated in 4.6.1 or 4.6.2.

The function center_coords is intended to shift the COM of the system to
bX/2, bY/2, 0 as described in the help documentation for bCenter. This is
accomplished by calculating a shift vector on line 153 (line numbers may be
off, sorry):

rvec_sub(box_center, com, shift);

and then subtracting the shift from all coordinates in a loop on line 162:

for (i = 0; (i < atoms->nr); i++)
   rvec_dec(x0[i], shift);

This is incorrect; the shift vector should be added to all coordinates, not
subtracted, for the COM to be centered since rvec_sub sets shift to
box_center - com. This can be fixed by changing rvec_dec to rvec_inc, which
I have verified returns the correctly centered coordinates.

The other issue is that the COM is computed by weighting by the mass of
each atom obtained from the topology, as shown in the command:

mm     = atoms->atom[i].m;

This is fine if the -dens flag is set to mass or electron; but for number
or charge the atom[i].m variable in the topology is reset in the function
g_density to either 1 (for -dens number) or the atom's charge (for -dens
charge). The system is thus centered not by the center of mass, but rather
by the geometric center/center of charge respectively. This could be fine
but may not be obvious to the user.

Again I added a report of both issues to the redmine report and wanted to
inform users of their existence in case this affects the use of the
g_density tool.

Best,
Reid