Analyzing MD from GROMACS

Here I do a short explanation on how to extract some magnitudes interesting for an analysis of you MD results:

Extract configurations
Extract movie
Extract Energy
Calculate Radial Distribution Function
Calculate dipole correlation
Angle Correlation Function (ACF)
Calculate S(q)

Goal: extract confs to be analyzed with ANGULA (gmx trjconv)

First of all, let's have a look on how to extract some confs in PDB format to be analyzed with ANGULA. This is quit easy: we will extract configurations starting at 100 ps and ending at 300 ps. Configurations will be saved each 5th frame.

gmx trjconv -f traj.trr -o system.pdb -pbc whole -b 100 -e 300 -skip 5 -sep

Goal: Movie with the trajectory (gmx trjconv)

Well, let's have a look, first of all, at the movie you got... it is nice. First of all let's translate the traj.trr file (where you have the trajectory of your molecules) to a PDB file, so that we can read it with VMD:

gmx trjconv -f traj.trr -o trajout.pdb -pbc whole

We have added -pbc whole so that molecules are kept together, and do not break because of periodic boundary conditions. So, get popcorn and soda and you should be able to watch the movie...

Goal: investigating energy, or any scalar (gmx energy)

First of all you have to extract the energy (or any other scalar that offers GROMACS to you). To do that, typing:

gmx energy -xvg none

is enough. It will generate an energy.xvg file Now we want to plot it. With the option -xvg the output format is simply x y (no titles) I have always used gnuplot, so if you do not specify -xvg none, the output format is for grace.

Goal: Calculating a Radial distribution Function (gmx rdf)

Well, well, well... radial distribution function. this is done in two steps:

You select the atoms you want to calculate the RDF from
You calculate the RDF

It could be easier, but it is not. To do that, imagine you have a molecule called MOL and an atom called AT (as wrote in the PDB or GRO file) You first do the index file using a configuration, for example the final one from a simulation (confout.gro):

gmx make_ndx -f confout.gro

Then a dialog will appear... am example will be more explanatory then the help of the command. If you want to select atom AT from molecule MOL you type

r MOL & a AT

... or maybe you are interested in all atoms that start with a C (all carbon atoms, for example). Then you would want to write

r MOL & a C*

This will generate an index group with the numbers of atoms AT in your molecule named MOL. The program will tell you the number of the newly created group (let's say 5). Then if you want to calculate the RDF from AT to AT, you just have to type

gmx rdf -f traj.trr -n index.ndx -ref 5 -sel 5 -o g_AT_AT.xvg

And that will generate an RDF called g_AT_AT.xvg

Goal: calculating dipole-dipole correlation function (gmx dipoles)

The goal is to calculate the dipole-dipole correlation function as a function of time: either that of the whole system or the average for all molecules.

[-corr total] In the first case we calculate [M(t)·M(0)]/M²(0) where M(t) is the dipole of the whole system.
[-corr mol] In the second case we calculate μ(t)·μ(0)/μ²(0) and do the average.

If you want to calculate the dipole correlation function for molecules (second case)

gmx dipoles -corr mol -P 1 -f traj.trr -c dipcorr.xvg

However this option is, of course, much more powerful than that. Look gromacs manual to learn more and extract the full power and potential of this option _{(this last sentence should be read as an energy drink commercial)}

Goal: Calculate the Angular Correlation Function

The question here is to take a vector, and then calculate the correlation of this vector, with itself... after some time t. More into specifics, the scalr product as in the previous case... but now for an arbitrary vector. You have two options:

Calculate the correlation of a vector coming from two points of the molecule. And therefore defined as i,j
Calculate the correlation of a vector coming from a plane. In this case three points i,j,k must be given!

To do that first we will create an index file with the points to be calculated, and then we will launch the Autocorrelation program.

STEP 1: selecting the atom pairs otr triplets:

In the first case you should be writing something like " a C1 | a C2" to generate an index file with atoms C1 and C2 for each molecule (i.e. atoms C1 or C2 i the selection ;-)
In the second cas (planes), it comes as no surprise, that you should be writing something as "a C1 | a C2 | a C3"

STEP 2: calculating the ACF

In the first case you must write something like "gmx rotacf -d". Note the -d switch: you are telling gmx that she will find only atom pairs. When you press enter, select the group with pair of atoms.
In the seconf case, the same without the switch -d...

Goal: calculating the total scattering function S(Q)

Dedicated to Neutron Lovers! : Follow this link to do that