[gmx-users] Re: nonormalize option in g_rotacf

[Andrew DeYoung]

Hi Justin,

Thanks so much.  I really appreciate your time looking into this.  Like you
said, it is hard to imagine that this could possibly be a bug.  Anyway, I
will submit an issue in a few days if there are no other comments.

Andrew

------------------
Date: Sun, 12 Jun 2011 20:43:39 -0400
From: "Justin A. Lemkul" <jalemkul at vt.edu>
Subject: Re: [gmx-users] nonormalize option in g_rotacf

Andrew DeYoung wrote:
> Hi, 
> 
> I am having some seeming difficulty with the -normalize option in
g_rotacf,
> which calculates rotational autocorrelation functions of molecules.  
> 
> My system consists of 254 SPC/E water molecules.  I equilibrated (3 ns, md
> integrator) and ran production dynamics (2 ns, md integrator) in the NVT
> ensemble, using Nose-Hoover temperature coupling.  I chose the cubic box
> size such that the density of the system is set to ~1 g/cm^3.  After I
> completed my production run, I used following command:
> 
> g_rotacf -f nvt-md.trr -s nvt-md.tpr -n allatoms.ndx -P 1 -b 0 -e 10 -o
> finalresults/normalized_10ps.xvg
> 
> where allatoms.ndx specifies all 762 atoms in triplets of atoms, with each
> triplet corresponding to one water molecule.  The default setting for
> -normalize is yes, so the above command should give me the _normalized_
> rotational autocorrelation function, averaged over all molecules.  I do
> indeed get a normalized ACF, which has C(t=0)=1 by definition.  (Recently
> (http://lists.gromacs.org/pipermail/gmx-users/2011-June/062005.html),
Justin
> kindly pointed out to me that, per the normalize_acf() function in
> src/tools/autocorr.c, the normalization is such that the ACF is 1 at time
> t=0.)
> 
> But, on the other hand, when I use the option -nonormalize, using the
> following command, I get a seemingly strange result:
> 
> g_rotacf -f nvt-md.trr -s nvt-md.tpr -n allatoms.ndx -nonormalize -P 1 -b
0
> -e 10 -o finalresults/unnormalized_10ps.xvg
> 
> When I execute this command, I obtain an ACF that is 1 at time t=0.  This
> cannot really be the correct unnormalized result, unless by some remote
> coincidence the unnormalized ACF happens to be 1 at time t=0.  
> 
> I have posted my results at
> http://www.andrew.cmu.edu/user/adeyoung/june12/june12.pdf .  If you have
> time, could you please help me think where I may have gone wrong?  At
first,
> I was convinced that I mixed up the files or did something else silly.
But,
> I've redone it, carefully naming, opening, and plotting the results, and I
> get the same: both the normalized and the unnormalized ACFs are 1 at time
> t=0, and the coordinates in the .xvg files are identical.  I am sure that
I
> am doing something wrong, but I have not found it yet.  Can you please
help
> me think what to try next?  I am using Gromacs version 4.5.4.
> 
> When I want to turn off normalization of the ACF, should I use
> "-nonormalize" (as I have done above, and as seems to be suggested by the
> manual), or should I use "-normalize no"?
> 

You're invoking the command correctly.

Your results indeed suggest that the normalization is being done regardless
of 
your command.  I looked a bit into the code and it seems that the
acf.bNormalize 
variable is initialized to TRUE in the function add_acf_pargs() in
autocorr.c; I 
see no later instance where this value can be changed to FALSE.  I did a
quick 
test where I changed the initial value of acf.bNormalize variable in
autocorr.c 
to FALSE and then ran g_rotacf with -normalize and -nonormalize.  In this 
instance, I got a more intuitive result - normalization in the case where it
was 
requested, and raw values with the same ratios as the normalized output in
the 
case where -nonormalize was used.

This very basic test suggests that normalization is not handled correctly,
but I 
will say that this is the first time I'm looking at the code (which is quite

complex), so it is very possible that I'm missing something.

I will also state the caveat that changing the value of acf.bNormalize is
just 
some quick and dirty idea that came to mind.  I do not suggest that you (or 
anyone else who may come across this thread!) do this and then base a whole
lot 
of results off of the resulting output without a developer saying otherwise.

The mechanism of command line parsing and ACF calculations are beyond what
I'm 
qualified to comment on authoritatively.  Hopefully someone who is
knowledgeable 
about such mechanisms will comment.  If not, please file an issue on 
redmine.gromacs.org about this unexpected outcome.  It smells like a bug,
but it 
seems that it should have affected all ACF results for all tools that
calculate 
ACF's.  Weird to have not heard about it until now, but I can't see any
other 
expectation based on your results.

-Justin

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

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

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