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" <jalem...@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 ======================================== -- gmx-users mailing list gmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
