[gmx-users] Rotational autocorrelation function of molecules using g_rotacf

[Andrew DeYoung]

Hi,

I am new to Gromacs and to the field of molecular dynamics, and I would like
to calculate the rotational autocorrelation function of molecules using
g_rotacf (http://manual.gromacs.org/online/g_rotacf.html).  

My system consists of 254 water molecules (so it's a very tiny system).
When I look in the command specification for g_rotacf in the manual, I see
that "three atoms (i,j,k) must be given in the index file, defining two
vectors ij and jk.  The rotational ACF is calculated as the autocorrelation
function of the vector n = ij x jk, i.e. the cross product of the two
vectors.  Since three atoms span a plane, the order of the three atoms does
not matter."

My molecules are water molecules, which of course each consist of three
atoms, hydrogens H1 and H2 and oxygen O.  My .gro structure file looks like
this:

Protein
  762
    1SOL     OW    1   1.798   1.179   0.524 -0.1733  0.7269  0.1573
    1SOL    HW1    2   1.754   1.268   0.520 -0.7587  0.4336  0.0304
    1SOL    HW2    3   1.822   1.149   0.432  0.0246  0.7483  0.2000
    2SOL     OW    4   0.736   0.930   1.977  0.3864 -0.9855 -0.7272
    2SOL    HW1    5   0.769   0.926   1.882  1.9432 -0.3254 -0.2249
    2SOL    HW2    6   0.814   0.938   2.039 -0.6975  0.1296  0.5182
    3SOL     OW    7   0.589   1.155   0.059  0.3331  0.4782  0.9113
    3SOL    HW1    8   0.611   1.206  -0.024 -3.4527 -0.6840 -0.9167
    3SOL    HW2    9   0.626   1.063   0.052 -0.1206  0.3293  0.4107
...

Of course, water is not protein!  This incorrect title is because my
original .pdb file incorrectly had the "Title" option set to "Protein".

Anyway, I decided to create an index file specifying atoms 2, 1, and 3 (HW1,
OW, and HW2), although the order does not matter.  In make_ndx, I used the
command "a 2 1 3" and saved it in an index file called water1.ndx.  Then I
used the following command to get the rotational ACF corresponding to the
first-order Legendre polynomial:

g_rotacf -f nvt-md.trr -s nvt-md.tpr -n water1.ndx -P 1 -nonormalize -o
rotacf/rotacf.xvg

and selected the new "a 2 1 3" option that I had previously created in the
index file.  I got a reasonable-looking plot.  I also calculated the
unnormalized version by omitting the -nonormalize option, and these two
plots turned out to be exactly identical, curiously.  

So, unless I have made some big error, I have calculated the rotational ACF
of water molecule 1 in my system, using a reference vector that is
perpendicular to the plane containing the three atoms.  But, looking briefly
in the literature, I see a lot of mentions of terms like "averaged
reorientational autocorrelation functions" and "averaged rotational
autocorrelation functions".  I am not certain, but perhaps this would mean
in my case that I should compute the rotational ACF of EACH water molecule
(using the respective reference vectors perpendicular to the atoms), and
then average the results, and finally plot the average results as a function
of time.  

It is not clear to me how I should do this.  The entry for g_rotacf in the
manual clearly states that three atoms must be specified in the index file.
I don't see how it can be possible to specify more than three atoms,
especially since four or more points do not in general constitute a plane.
But just out of curiosity, I made a new index file, this time specifying the
atoms of the first two water molecules, using "a 2 1 3 5 4 6".  I then ran
this command: 

g_rotacf -f nvt-md.trr -s nvt-md.tpr -n twoindex.ndx -P 1 -nonormalize -o
rotacf/tworotacf_unnormalized.xvg 

I expected to get an error message or a completely nonsensical result (or
the identical result as before, where the "extra" atoms 5 4 6 ignored by the
algorithm), but, suprisingly to me, I got a result similar to, but not
identical to, my previous "one molecule" calculation.  I have posted these
results at http://www.andrew.cmu.edu/user/adeyoung/june7/rotationalacf.pdf
in case it is helpful, where the "one molecule" case is in red and the "two
molecule" case is in blue.

If you have time, can you please help me find more documentation or
explanation of g_rotacf or about rotational autocorrelation functions in
general?  I feel completely silly, because although g_rotacf is mentioned on
page 206 and specified on page 300 of version 4.5.4 of the manual, I feel
like I am just not seeing some of these points; I am sure they are there, I
am just not seeing it.

I am mainly looking for answers to these questions: 

(1) When I specify six atoms (2 1 3 5 4 6) instead of three (2 1 3), what am
I doing or what is g_rotacf doing with the "extra" atoms?  From what I have
read so far, I have absolutely no reason to believe that g_rotacf is
calculating the rotational ACF of atom sets 2 1 3 and 5 4 6 and then
computing the average.  

(2) Is there a way to calculate the average rotational autocorrelation
function?  Should I be thinking about writing a program to find separately
the rotational ACF of each of my 254 water molecules, and then average them
at the end?  

(3) What is the normalization being done by g_rotacf (when -nonormalize is
omitted)?  Is it just that every element in the rotational ACF list is
divided by the rotational ACF at time t=0, so that the initial value of the
rotational ACF is always 1?

Thank you very much for all your time!

Andrew DeYoung
Carnegie Mellon University