[gmx-users] Computing vectors normal (perpendicular) to molecules

[Andrew DeYoung]

Hi,

Suppose that I have a system of benzene molecules (in reality, my system is
more complicated, but for my question, it will be simpler to consider just
an ensemble of benzene molecules).  

I would like to find a vector normal (i.e., perpendicular) to the plane of
each benzene molecule in my system.  I want to track these "normal vectors"
over time; eventually, I would like to calculate a histogram of angles theta
between the vectors normal to the benzene molecules and some fixed
"laboratory" axis.

If you have time, can you please tell me if there is any Gromacs utility
that already does this -- i.e., computes normal vectors like this?  It does
not seem so, even with the functionalities provided by g_angle.

This may not be the most rigorous idea, but one way to do this might be to
compute the vector cross product of two vectors.  Suppose that a benzene
molecule has carbon atoms named C1, C2, C3, C4, C5, and C6 (arranged in that
order around the ring). Then I could, for example, find two vectors:

\vec{C2-C1}
\vec{C2-C3}

where \vec{C2-C1} is the vector from C2 to C1, and \vec{C2-C3} is the vector
from C2 to C3. 

The cross product of those two vectors gives a vector that is perpendicular
to both vectors.

Is there any way to compute cross products in the Gromacs utilities?  It
does not seem so, but I am not certain.

So, I am thinking that I will need to use g_traj to extract the x, y, and z
coordinates of the atoms of interest (i.e., C1, C2, and C3 in each of the
many benzene molecules in the system).  Then I will need to write a C or
Fortran script to find the relevant vectors and compute their cross
products, at every timestep in the trajectory.

However, to complicate things even more, I may in the future wish to make
dynamic selections of benzene molecules, using g_select.  I may want to
consider only the benzene molecules satisfying z<12, for example.  Since the
benzene molecules obviously move, the particular benzene molecules that I
consider for the vector/cross product calculation will change over time. I
am able (using the -oi option in g_select in conjunction with a script of my
own) to generate an index file with the indices (whose groups correspond to
the simulation timestep number) of the C1, C2, and C3 atoms in benzene
molecules whose centers-of-mass satisfy z<12.  But is there any way I can
feed this index file to g_traj without having to call g_traj as many times
as I have timesteps in my trajectory?

Thanks so much for your time!

Andrew DeYoung
Carnegie Mellon University