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

[Justin A. Lemkul]

Andrew DeYoung wrote:
> 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.  
> 

That's what it's doing.  The documentation should indicate that the index groups 
should contain atom triplets, rather than simply saying "three atoms."  I'll fix 
this.  You can see all the inner mechanics in src/tools/gmx_rotacf.c.

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

Not necessary.  This is what g_rotacf is doing.  Note that this is (very 
obscurely) implied by the default setting of -aver (which is "yes," meaning it 
will average over all molecules).

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

Yep, per the normalize_acf() function in src/tools/autocorr.c.

-Justin

> Thank you very much for all your time!
> 
> Andrew DeYoung
> Carnegie Mellon University
> 

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================