[gmx-users] Important bugfix for g_hbond

[Erik Marklund]

Hi,

It's a matter of how you count the hbonds. One water can participate in 
four hydrogen bonds (sort of); two as an acceptor, and two as a donor. 
But that leads to double counting of hbonds if we think that Nmax = 
Nwater*4. Analogously: if a long line of people are connected by holding 
hands, every person make two connections, but the total number of 
connections divided by the number of poeople approaches one as the 
number of people approaches infinity. Note how in the case of your hbond 
calculations 1.8*2 = 3.6.

Cheers,

Erik

Baofu Qiao skrev 2011-05-17 15.57:
> Hi Erik and others,
>
> After the bugfixed, I still have one problem: the average HB per SPC/E 
> water is calculated to be ~1.8 in a pure water system from g_hbond, 
> which has been previously reported to be 3.6 (Kumar, Schmidt & Skinner 
> JCP, 2007, 126, 204107) . I have no idea where the difference comes 
> from. Could you provide some help? Thanks a lot!
>
> Some info:
> ******************************************
> 4000 SPC/E water simulated at T=298K using NTP;
> "Average number of hbonds per timeframe 7217.386 out of 8e+06 
> possible" was reported by (echo -e "SOL\nSOL"| g_hbond2 -f T298 -s 
> T298 -n -merge). Note that the value of 7217 is independent of 
> -merge/-nomerge.
> I am using Gromacs 4.5.3 with your 'bugfix". Gromacs 4.0.7 give the 
> similar value of about 7210 with -merge or -nomerge.
> ******************************************
>
> Any suggestion is appreciate!
>
> Best wishes,
> Baofu Qiao
>
> On 05/10/2011 03:08 PM, Erik Marklund wrote:
>> Hi,
>>
>> There have been reports about inconsistencies between older (<= 
>> 4.0.7?) and newer versions of g_hbond, where the older seem to have 
>> been more reliable. I found and killed the bug that caused the newer 
>> versions to miscount the hbonds. Checkout release-4-5-patches to get 
>> the bugfix, or patch it yourself by commenting out line 1497, which 
>> reads "return hbNo;":
>>
>>      if (bBox){
>>         if (d>a && bMerge && (bContact || isInterchangable(hb, d, a, 
>> grpd, grpa))) { /* acceptor is also a donor a\
>> nd vice versa? */
>>             /* return hbNo; */
>>             daSwap = TRUE; /* If so, then their history should be 
>> filed with donor and acceptor swapped. */
>>         }
>>
>> Simple as that.
>>
>
>


-- 
-----------------------------------------------
Erik Marklund, PhD student
Dept. of Cell and Molecular Biology, Uppsala University.
Husargatan 3, Box 596,    75124 Uppsala, Sweden
phone:    +46 18 471 4537        fax: +46 18 511 755
erikm at xray.bmc.uu.se    http://folding.bmc.uu.se/