Dear all,
This is a follow-up to an old thread on g_principal, which continues (as of
version 4.6.5) to suffer from what I would call a bug or, at least, a very
bad
and misleading choice of output file names. This message is primarily a
further
warning to users, but I also hope that it promotes the solution of the
problem
(I explicitly indicate below the small required code changes).
The program g_principal "calculates the three principal axes of inertia for
a
group of atoms" and generates three output files with the default names
axis{1,2,3}.dat. Although the content of these files is not explained,
their
names naturally suggest that axisN.dat contains the xyz components of the
Nth
principal axis (presumably in the conventional major-to-minor order).
Indeed,
each of these files contains (besides the time) 3 real values per frame,
and
this interpretation also yields three vectors which are orthonormal, as one
would expect for the new frame defining the principal axes.
However, this "natural" interpretation is wrong, and not because of a
different
axes order. As correctly identified by Chris Neale in the message included
below, the file axis1.dat contains the x components of the 1st, 2nd and 3rd
axes, the file axis2.dat their y components, and the file axis3.dat their z
components (he checked with VMD, we did it with in-house programs in C and
Octave). I think this a rather convoluted and extremely misleading way to
output
the principal axes.
My guess is that this problem derives from a simple confusion between the
form
of the orthogonal matrix obtained when solving the eigenvalue problem --
depending on the adopted algorithm, the vectors defining the new basis (the
principal axes) can be either the columns or the rows of this matrix. I
didn't
look beyond gmx_principal.c, but the problem can be easily solved there by
replacing the current lines
fprintf(axis1, "%15.10f %15.10f %15.10f %15.10f\n", t,
axes[XX][XX], axes[YY][XX], axes[ZZ][XX]);
fprintf(axis2, "%15.10f %15.10f %15.10f %15.10f\n", t,
axes[XX][YY], axes[YY][YY], axes[ZZ][YY]);
fprintf(axis3, "%15.10f %15.10f %15.10f %15.10f\n", t,
axes[XX][ZZ], axes[YY][ZZ], axes[ZZ][ZZ]);
with
fprintf(axis1, "%15.10f %15.10f %15.10f %15.10f\n", t,
axes[XX][XX], axes[XX][YY], axes[XX][ZZ]);
fprintf(axis2, "%15.10f %15.10f %15.10f %15.10f\n", t,
axes[YY][XX], axes[YY][YY], axes[YY][ZZ]);
fprintf(axis3, "%15.10f %15.10f %15.10f %15.10f\n", t,
axes[ZZ][XX], axes[ZZ][YY], axes[ZZ][ZZ]);
With this change, the file names would make perfect sense, with axisN.dat
simply
containing the components of the Nth principal axis. (Note that both the
by-row
and by-column readings give orthonormal vectors, because this is an
orthogonal
matrix.)
The file moi.dat also produced by g_principal is fine, containing the
moments of
inertia along the principal axes in the proper order (lowest to highest,
since
the inertia _around_ those axes increases as one goes from the major to the
minor).
I believe that, as it stands now, g_principal is misleading many users into
the
wrong interpretation of its output. Maybe some developer wants to have a
look at
this issue and introduce my suggested fix.
