[gmx-users] Average no of hbonds

[Justin A. Lemkul]

Steven Neumann wrote:
> 
> 
> On Tue, Sep 6, 2011 at 3:09 PM, Justin A. Lemkul <jalemkul at vt.edu 
> <mailto:jalemkul at vt.edu>> wrote:
> 
> 
> 
>     Steven Neumann wrote:
> 
> 
> 
>         On Tue, Sep 6, 2011 at 2:50 PM, Mark Abraham
>         <Mark.Abraham at anu.edu.au <mailto:Mark.Abraham at anu.edu.au>
>         <mailto:Mark.Abraham at anu.edu.__au
>         <mailto:Mark.Abraham at anu.edu.au>>> wrote:
> 
>            On 6/09/2011 11:37 PM, Steven Neumann wrote:
> 
>                Dear Gromacs Users,
>                I am calculating hbonds between my 10 ligands and each
>                residue... How does Gromacs calculate average number of
>         hbonds
>                per timeframe?
>                Example:
>                for Glycine:
>                Av. num of hbonds/timeframe
>                0.96
>                To check whether that is correct I added all hbonds
>         formed with
>                Glycine during my simulation time over 2000 (each 50ps)
>                timeframes and I obtained value: 2900. If you divide it
>         by 2000
>                you will never get 0.96 obviously. How does Gromacs
>         calculate it?
> 
> 
>            Presumably by adding up the number in each time frame and
>         dividing
>            by the number. Unless you have a highly-exposed glycine, forming
>            more than one H-bond seems unlikely, and your 2900 number
>         suggests
>            it's happening lots of times. I expect you're comparing a sum of
>            oranges with a sum of apples, but without more information about
>            what you're doing in your attempt to check, we can't help much.
>            Please be sure to read g_hbond -h, and the legends of any
>         .xvg files
>            you're looking at.
> 
>            Mark
>            Thank you Justin and Mark,
> 
>          Yes, just found the mistake... I added pairs within 3.5 A
>         instead of hbonds so the Gromacs calculation is correct. So the
>         value is 0.96 for Glycine. So in this case it is a strong
>         interaction. Are there any criteria to assess strenght of
>         interactions or binding affinity for residues? Was it described
>         anywhere for each residue to specify these values with respect
>         to possible hbonds for each residue (e.g. Glycine - 3)?
>          
> 
> 
>     Binding affinity can be determined from free energy calculations,
>     but for the case of multiple ligands this would be an incredibly
>     complex calculation (or series of calculations, really).  You might
>     be able to make some argument about occupancy of available hydrogen
>     bonding sites.
> 
>     -Justin
> 
>  
> So will it be possible to exteract one ligand which from my 
> visualisation occupy one residue with highest number of av. hbonds 
> (others aggregate on bounded or bind weaker) and then calculate free energy?
>  

Something like that might be possible, by calculating free energies over an MD 
rerun, but I've never tried it so I won't speculate.  Consult the literature for 
precedent and procedure and have a look at the free energy tutorial on the 
Gromacs site to see if you can piece it all together.  It won't be a trivial 
exercise.

-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

========================================