Bug #417

ACF incorrectly calculated w Luzar/Chandler model

Added by Erik Marklund over 9 years ago. Updated about 9 years ago.

analysis tools
Target version:
Affected version - extra info:
Affected version:


The ACF is calculated for every existence function minus the average number of hbonds, i.e. C(t) = <(h(0)-<h>) * (h(t)-<h>)> / N, where N is for normalization. I believe this is wrong. In Luzar J. Chem. Phys. 113 (2008) the average is not subtracted like that. Subtracting <h> totally ruins the probability interpretation of the ACF.

I'll have a go on this bug myself.


#1 Updated by Erik Lindahl over 9 years ago

Has this been solved yet?

#2 Updated by David van der Spoel over 9 years ago

This was a patch for not treating periodicity correctly. If we treat periodicity as we should (do we?) subtracting the <h> can be kicked out.

#3 Updated by Erik Marklund over 9 years ago

I disagree. I really think this change should be done regardless of how (non-)periodicity is treated. I think I did the change, but didn't have time to test it thoroughly before the summer holidays.

I just got back to civilized areas with proper internet connections, so expect better communication from my part from now on. I'll be back at the lab some time next week and resolve this.

David, could you please guid me to where <h> is subtracted in Luzar and Chandler's papers? As I see it it goes against their equations.

#4 Updated by David van der Spoel over 9 years ago

It is not. I made it up.

#5 Updated by Erik Lindahl over 9 years ago

Any update on this?

If we don't have any obvious fix I would strongly advocate that we should follow the equations in the paper we cite by default. If somebody has a particular case where they need to subtract the average, can't that be done manually in post-processing?

#6 Updated by David van der Spoel over 9 years ago

I discussed it with Erik M. We agree that it should be thrown out. Erik will try to do it (but he's doing a course right now).

#7 Updated by Erik Lindahl about 9 years ago

I'm sorry, but the code in gmx_hbond() is simply way too complex for me to fix this without extensive testing and debugging, and several parts of it isn't even supported without GSL.

Is this simply the subtraction that present occurs on line 2811? It seems somebody familiar with the code should be able to solve this in 30 seconds?

Unfortunately I think this another reason it is time to move all the advanced functionality into a contrib tool, and recreate a much simpler g_hbond from scratch, which only does the primary data extraction and leaves the analysis to the user.

#8 Updated by David van der Spoel about 9 years ago

I fixed it anyway in commit 24f544ad127047b6dd435e70c2b81700c89fb72c.
Further development after 4.5.

#9 Updated by Erik Lindahl about 9 years ago

Great, thanks! I didn't dare to, since I wasn't 100% sure. In that case we can close it completely for now.

Also available in: Atom PDF