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

Antonio Baptista baptista at itqb.unl.pt
Sat Feb 8 00:18:38 CET 2014

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]);


         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.


Antonio M. Baptista
Instituto de Tecnologia Quimica e Biologica, Universidade Nova de 
Av. da Republica - EAN, 2780-157 Oeiras, Portugal
phone: +351-214469619         email: baptista at itqb.unl.pt
fax:   +351-214411277         WWW:   http://www.itqb.unl.pt/~baptista

Nov 02, 2010; 4:46pm Chris Neale
Principal axis (g_principal)

Dear Sarah,

I just ran into the same thing so I'm posting this as a very late
response to your querry.

I think that g_principal prints the transpose of the principal axes
matrix, which is the exact opposite of what one would expect.

(1) First, I define a system that is longest in Z, second longest in
Y, and shortest in X:

cneale at iram:~$ cat a.gro
      1CL-     CL    1   0.000   0.000   0.000
      2CL-     CL    2   0.000   1.000   0.000
      3CL-     CL    3   0.000   2.000   0.000
      4CL-     CL    4   1.000   0.000   0.000
      5CL-     CL    5   1.000   1.000   0.000
      6CL-     CL    6   1.000   2.000   0.000
      7CL-     CL    7   0.000   0.000  10.000
      8CL-     CL    8   0.000   1.000  10.000
      9CL-     CL    9   0.000   2.000  10.000
     10CL-     CL   10   1.000   0.000  10.000
     11CL-     CL   11   1.000   1.000  10.000
     12CL-     CL   12   1.000   2.000  10.000
50 50 50
cneale at iram:~$ cat a.top
#include "ffoplsaa.itp"
#include "ions.itp"

[ system ]
; name

[ molecules ]
; name        number
CL-  12


(2) Next I calculate the principal axes:

$ echo 0 |g_principal -f a.gro -s a.tpr >/dev/null 2>&1; cat axis1.dat
axis2.dat axis3.dat
     0.0000000000        0.0000000000     0.0000000000     1.0000000000
     0.0000000000        0.0000000000     1.0000000000     0.0000000000
     0.0000000000        1.0000000000     0.0000000000     0.0000000000


(3) Finally, I rotate it by 90 deg and calculate the principal axes:

$ editconf -f a.gro -o b.gro -rotate 90 0 0 >/dev/null 2>&1; echo 0
|g_principal -f b.gro -s a.tpr >/dev/null 2>&1; cat axis1.dat
axis2.dat axis3.dat
     0.0000000000        0.0000000000     0.0000000000     1.0000000000
     0.0000000000        1.0000000000     0.0000000000     0.0000000000
     0.0000000000        0.0000000000     1.0000000000     0.0000000000


(4) I use VMD to confirm that the longest axis is in Y. Thus one must
read the columns in axis1.dat as time1, X1, X2, X3 whereas one would
intuit that they read as time1, X1, Y1, Z1.

I didn't test this extensively, but if my colleague determines that I
am incorrect in what I have posted here then I will post that to this

If I am correct, then it would be very useful to get some #titles in
the axisN.dat output files from g_principal.



Dear gmx users,

I am doing several analyses (version 4.0.4) on my simulations with 
organic molecules inserting into a DMPC bilayer. Now I would like to 
whether the small molecule inserts into the membrane with a specific 
angle to
the membrane normal (or z-axis). I have used two approaches:

1) My molecule is a cyclohexene ring with two "para" substituents. I 
made an
index group with carbon number 3 and number 10 (in different groups) 
corresponds to the first substituent-atoms on each side of the ring. 
vector describes the shape of the molecule quite well. I then used 
g_bundle to
calculate the angle between this vector and the z-axis by echoing first 
C3-group and then the C10-group:

echo 0 1 | g_bundle -f xxx.xtc -s xxx.tpr -n xxx-vector.ndx -na 1 -z 

This gives fairly reasonable results and angle distribution.

2) I also calculated the principal axis by g_principal where I echoed 
small molecule as a whole:

echo 0 | g_principal -f xxx.xtc -s xxx-1.tpr -n xxx-lim.ndx -a1
xxx-principal.dat -a2 xxx-axis2.dat -a3 xxx-axis3.dat

To get the angle between the principal axis and the z-axis one takes 
the arcos
the the z-component. This created quite a different and very broad 

This puzzled me very much since the difference between the C3-C10 
vector and
the principal axis should not be that great.

I then used a small python program to calculate the principal axis of 
molecule and the result was again quite different than from
g_principal. I have
tried to rotate each of my substituent groups by 120 degrees (creating 
different conformations) and then calculate the difference in the 
axis but this creates only a deviation of max. 16-17 degrees not at all
accounting for the very different angle distributions seen. I have not
numbers but I can do so along with the python code. I would like to ask 
anyone can see any obvious mistakes or know of something special to be
aware of
when using g_principal?

Thank you,


More information about the gromacs.org_gmx-users mailing list