[gmx-users] g_principal -- bug or very bad choice of filenames

[Justin Lemkul]

On 2/7/14, 6:08 PM, Antonio Baptista wrote:
> 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.
>

Please file an issue on redmine.gromacs.org with all of the above information. 
Thanks for the thorough report!

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 601
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==================================================