[gmx-users] How Gromacs h_bond module compute hydrogen correlation function ?

[Erik Marklund]

Hi again,

I had a look at your post at ResearchGate, where you ask the same question and provide both source code and output from your. Two things differ between gromacs’ calculation and yours. First of all, gromacs normalises the correlation function so that it is 1 at t=0. This is common practice and allows for a probability interpretation of the ACF. Secondly, you may have noticed that the gromacs output contains several datasets, one of which is the ACF with a correction that compensates for finite size effects, and one "raw” ACF that resembles yours, apart from said normalisation. Note that the ACF will not approach zero in a periodic system, since pairs that would diffuse apart in a large system will easily find each other again because the diffusion wraps around the pbc. Hope that helps.

Kind regards,
Erik

> On 14 Mar 2016, at 10:22, Erik Marklund <erik.marklund at chem.ox.ac.uk> wrote:
> 
> Hi,
> 
> Gromacs uses FFTs to calculate C(t), exploiting that convolutions (such as an autocorrelation) turns into simple muliplications in Fourier space. If you are interested in the details, have a look at gmx_hbond.c. In the function do_hbac(), look for the last call to low_do_autocorr() and the code around it. It should be under “case AC_LUZAR:”. There may be other calls to low_do_autocorr(), depending on your version, but they concern other kinetic models and alternative bond definitions that aren’t fully supported.
> 
> Kind regards,
> Erik
> 
>> On 12 Mar 2016, at 07:13, 백호용 (자연과학부) <ghdyd158 at unist.ac.kr> wrote:
>> 
>> Dear gromacs users and developers, I have a question about the methodology that gromacs h_bond module compute H_bond correlation function, C(t) from trajectory.
>> 
>> I want to compute hydrogen bond life time between a carbonyl oxygen of single Etoac and water. For example, I made a system which contain a single Etoac molecules and about 2000 water molecules. then I used h_bond module to get hbmap.xpm and hbond.ndx and I analyzed those two so that I can get existence function of each carbonyl oxygen - water pair. In detail, during 500 ps simulation, 23 water molecules formed hydrogen bond with carbonyl oxygen of Etoac at least one time. So, I computed 23 existence functions, each having 500 ps length. Then I use those 23 set of existence function to compute h_bond correlation function.
>> 
>> I'll upload Matlab code (corr_031116.m )with which I computed correlation function and please tell me what the problem is..
>> 
>> So, What I want to know is that How Gromacs practically compute Correlation function, in detail..
>> 
>> Thank you for reading my email and Thank you for your reply in advance..
>> 
>> -- 
>> Gromacs Users mailing list
>> 
>> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>> 
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>> 
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.
> 
> -- 
> Gromacs Users mailing list
> 
> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
> 
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> 
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.