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

[Mark Abraham]

On 2/07/2012 6:50 AM, Andrew DeYoung wrote:
> 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.

Yes there are tools that can do this, e.g. g_sgangle or perhaps others. 
Manual 7.4 is a good place to start your search, then to appendix D.

> 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?

Not until GROMACS 5.

Mark