SV: SV: SV: [gmx-users] Hydrogen bonding

[Justin A. Lemkul]

Sarah Witzke wrote:
> 
> 
> -----Original Message----- From: gmx-users-bounces at gromacs.org on behalf of
> Justin A. Lemkul Sent: Thu 19-11-2009 02:09 To: Gromacs Users' List Subject:
> Re: SV: SV: SV: [gmx-users] Hydrogen bonding
> 
> 
> Sarah Witzke wrote:
>> 
>> 
>> Sarah Witzke wrote:
>> 
>> <snip>
>> 
>>> ARRGH, I'm sorry, things went too quick :-( So the lifetime in average
>>> per hbond is 628.571 ps?
>>> 
>> Yes, per the calculation.  For a bit more about the analysis, see the
>> "Please read and cite" notices, as well as this thread:
>> 
>> http://lists.gromacs.org/pipermail/gmx-users/2009-February/039955.html
>> 
>>> The note about the correlation function would imply that the correlation 
>>> function itself has not converged until its value is < 0.001.  This is
>>> usually a result of insufficient data, either the length of the
>>> simulation, or number of frames analyzed (based on the spacing of the
>>> frames).
>>> 
>>> Hmm, so using 120 ns simulation from a standard .xtc file (timestep of 10
>>> ps) is not long enough. Perhabs I should try the original .trr file...
>>> The .xtc file I'm using only contains DMPC and the small molecule, no
>>> solvent - I can't see this would make a difference in this case, am I
>>> right?
>>> 
>> The correlation will depend on how much the interactions are changing over
>> the period you analyzed.  If you are analyzing a small molecule and DMPC,
>> water should not matter.
>> 
>> I have now tried with the full length (220 ns) .trr file and the -ac output
>> still prints a warning:
>> 
>> ACF 106/106 Normalization for c(t) = 1.23261 for gh(t) = 5.58815e-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.000172327 +/- 0.00277407
>> 
>> In your opinion, does this mean that I cannot trust the value of the
>> lifetime, since the correlation function is not converging?
>> 
> 
> I'm not entirely sure.  Please see the posts from David in the thread I 
> referenced before about some potential issues with the ACF in the g_hbond 
> calculation.  I assume you are using version 4.0.x?
> 
> I'm using version 4.0.4 for consistency. I have read the posts from you and
> David van der Spoel
> (http://lists.gromacs.org/pipermail/gmx-users/2009-February/039955.html) as
> well as the article given by g_hbond. I guess this auto correlation function
> and the lifetime derived by this is... well, I'm not sure I trust it enough
> to put it in a paper - especially since my acf doesn't converge (which I
> still find strange, it should not matter that I have only saved coordinates
> every 10 ps in my .trr file when using the Luzar acfs-approached, right?).
> 
> The - hbm option i g_hbond creates a matrix that can be converted to a .eps
> file by usin xpm2ps. The picture obtained is nice looking but lacking in some
> of the information I would like to derive. I would like to know how many
> frames each hbond exists (another way to consider lifetime). I have opened
> the .xpm file in a text editor and there are a long list of lines starting
> similar to tgis one: /* x-axis:  100000 100010 100020 Can I somehow obtain
> the number of frames for each hbond from this?
> 

Either write a script that counts the number of instances of "o" characters per 
line in the .xpm file, then divide by the total number of frames, or iteratively 
call g_hbond for each of the H-bonds that were initially identified (and whose 
indices are now written to hbond.ndx from -hbn).  The latter approach will still 
require some post-processing, of course, since each output .xvg file will be a 
series of either 0 or 1, indicating the absence or presence of that particular 
H-bond.

-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

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