[gmx-users] Question about documentation of -a switch in g_hbond

[Mark Abraham]

On 11/06/2012 12:41 AM, Justin A. Lemkul wrote:
>
>
> On 6/10/12 10:36 AM, Andrew DeYoung wrote:
>> Greetings,
>>
>> This might end up being a silly/embarrassing question, and if so, I
>> apologize.  I feel like I may be making a conceptual mistake, but I'm 
>> not
>> sure.
>>
>> Is it true that a hydrogen bond is of the following form?
>>
>> Donor---Hydrogen ... Acceptor
>>
>> Is this the correct order?  I think so.  For example, in my system, I 
>> have:
>>
>> Oxygen---Hydrogen ... Fluorine
>>
>> In the documentation, the -a option has this description:
>>
>> "Cutoff angle (degrees, Acceptor - Donor - Hydrogen)"
>>
>> Why is the order Acceptor - Donor - Hydrogen (or, equivalently, 
>> Hydrogen -
>> Donor - Acceptor)?  Why isn't "Hydrogen" instead between "Donor" and
>> "Acceptor"?  Which angle does -a specify?
>>
>
> Because it's easy to define from an algorithmic standpoint.  A 
> hydrogen bond is stronger depending on how close to linearity it is.  
> Measuring an acute angle given by -a says that if the A-D-H angle is 
> less than 35 degrees, you can vary between 0 and 35 easily to show the 
> same thing as being within a certain range of being 180 degrees in 
> terms of D-H-A.

The component of H-bond detection based on the cut-off angle is not 
implemented as a cut-off in angle space, but rather in the space of the 
cosine of the angle. Cosines are easily computed using the inner product 
of the vectors that form them. There's no real advantage to computing 
cos(A-D-H) vs cos(D-H-A) because in either case you just compare the 
value in the appropriate sense with the appropriate value of the 
(signed) cos(acutoff). Defining the H-bond in terms of A-D-H has the 
advantage that if you later want to consider the angle distribution, the 
relevant (signed) angle range is centered around zero degrees. D-H-A 
would have a cluster below 180 and a cluster above -180, which is less 
convenient to deal with.

For a given r(D-H), you can define the set of legal H-bond positions 
equivalently from either A-D-H or D-H-A, of course. Since r(D-H) is 
often constant for MD simulations, the conceptual difference between the 
definitions becomes moot.

In neither case do you want to compute the angle with acos later, 
because the numerical accuracy of acos is poor around 0 and +/-180 from 
its near-linearity. Unfortunately g_hbond does do this, so I think I'll 
patch it.

Mark