SV: [gmx-users] Hydrogen bonding

[Sarah Witzke]

Sarah Witzke wrote:
> Dear gmx-users,
>
>
>
> I have done simulations of one small molecule that diffuses into a DMPC
> membrane. This small molecule contains an alcohol group and is therefore
> capable of hydrogen bonding to the oxygens of DMPC (phosphate and glycerol
> region).
>
> I have read the manual (section 8.12 and g_hbond -h), searched the mailing
> list and google but I have not been able to find a more thorough description
> of the output possibilities than in the manual.
>
> I have tried three different approaches:
>
> 1. The -OH group of the small molecule and the glycerol oxygens
>
> 2. The -OH group of the small molecule and the phosphate oxygens
>
> 3. The small molecule and DMPC (no subgroups)
>
> No. 1 gives 38 hbond, no. 2 gives 15 hbonds and no. 3 gives 53 hbonds. So 1 +
> 2 = 3, which is fine.
>

I'm assuming this is just some theoretical limit that you have established, and
not something that has actually been calculated, correct?

This is from the output of the program, you have explained it below :-)

>
>
> Below is the output from no. 1: (gromacs 4.0.4)
>
> Specify 2 groups to analyze:
>
> Selected 0: 'O11_&_PALC_H12_&_PALC'
>
> Selected 1: 'O7_&_DMPC_&_DMPC_O9_&_DMPC_O10_&_DMPC'
>
> Checking for overlap in atoms between O11_&_PALC_H12_&_PALC and
> O7_&_DMPC_&_DMPC_O9_&_DMPC_O10_&_DMPC
>
> Calculating hydrogen bonds between O11_&_PALC_H12_&_PALC (2 atoms) and
> O7_&_DMPC_&_DMPC_O9_&_DMPC_O10_&_DMPC (384 atoms)
>
> Found 1 donors and 385 acceptors
>
> Making hbmap structure...done.
>
>
>
> Will do grid-seach on 15x15x24 grid, rcut=0.35
>
> Found 15 different hydrogen bonds in trajectory
>
> Found 23 different atom-pairs within hydrogen bonding distance
>
> Merging hbonds with Acceptor and Donor swapped
>
> - Reduced number of hbonds from 15 to 15
>
> - Reduced number of distances from 23 to 23
>
> Average number of hbonds per timeframe 0.083 out of 192.5 possible
>
>
>
> What does these "15 different hydrogen bond in trajectory" mean? I don't
> understand this. I also don't understand "Average number of hbonds per

That means quite literally what it says: there are 15 distinct hydrogen bonds
that form at some point in your trajectory.  They are listed in hbond.ndx
(output of -hbn) and mapped in hbmap.xpm (from -hbm).  You should find that
there are 15 :)
 
Yes, I see this. Do you by the way know why this .ndx files under the title [ donors_hydrogens_DMPC ] lists a lot of  non-heteroatoms (carbon atoms)?

> timeframe 0.083 out of 192.5 possible" - 192.5 possible hbonds?? Can anyone
> shed some light on this?

You have 385 H-bond acceptors (see the output above).  It appears that g_hbond
makes a simple assumption that if half of these were occupied (with the other
half of the group serving as a donor, which may or may not ever have any
physical significance), then the maximum number of H-bonds would be 385/2 = 192.5.

The average printed is simply the average number of H-bonds that are present
between your two groups at any given time.  If your small molecule has only one
-OH group that you're considering, the maximum number of H-bonds per timeframe is 1.
 
Ok.

>
>
>
> Another question relates to the lifetime of the hbond calculated when the
> "-life" flag is given. The produced .xvg file contain three columns: time,
> p(t), and t p(t). What is p(t) and t p(t)? And how can I find the lifetime?
>
>

The explanation is probably in one of the "Please read and cite the following
reference" messages that g_hbond spits out.  Looks like some sort of probability
function.  You should also be cautioned by this warning from using the -life option:

"Note that the lifetime obtained in this manner is close to useless
Use the -ac option instead and check the Forward lifetime"
 
When I use the -ac option I get this output: 
...
  ACF 53/53
  Normalization for c(t) = 1.19876 for gh(t) = 9.96115e-05 
  WARNING: Correlation function is probably not long enough
  because the standard deviation in the tail of C(t) > 0.001
  Tail value (average C(t) over second half of acf): 0.0155651 +/- 0.014574
  Hydrogen bond thermodynamics at T = 298.15 K
  Fitting parameters chi^2 = 5.76828e-05
  Q =          0
  --------------------------------------------------
  Type      Rate (1/ps) Time (ps)  DG (kJ/mol)  Chi^2
  Forward         0.001    682.571      20.705  5.76828e-05
  Backward        0.013     75.432      15.245
  One-way         0.001   1001.338      21.655
  Integral        0.000   3419.839      24.700
  Relaxation      0.001   1216.745      22.138
  100%
...
 
When looking at the produced .xvg file I see these three data sets as a function of time:
  @ s0 legend "Ac\sfin sys\v{}\z{}(t)"
  @ s1 legend "Ac(t)"
  @ s2 legend "Cc\scontact,hb\v{}\z{}(t)"
  @ s3 legend "-dAc\sfs\v{}\z{}/dt"

 
In neither the output nor the .xvg file I see any mention of 'Forward lifetime'. Perhabs this is because of the warning 'Correlation function is probably not long enough because the standard deviation in the tail of C(t) > 0.001' which I sadly don't understand. Is this a matter of actual length of the .xtc file? I use the last 120 ns of a 220 ns simulation. Any suggestions?
 
Thank you very much for all your help!

-Justin

--
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================
--
gmx-users mailing list    gmx-users at gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php