ACF function hbond - few questions

GROMACS version: 2018.4
GROMACS modification: No
I have some questions. I want to compute average hbond lifetime between lipids and water molecules. To do this I use ACF
For example gmx hbond -f eq3_900_1000.xtc -s eq3.tpr -b 900000 -e 950000 -ac eq3_ac.xvg

  1. I get time (0-5000 ps) and 4 values, but I can’t find good description. In my xvg file I have
    @ 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”
    This is my images (x - time from 0 to 5000 ps, and 4 different values)

    2.Why in the first value (“Ac\sfin sys\v{}\z{}(t)”) at the end I have negative number like -0.0175003
    The first two values are autocorrelation function, but third and fourth?

Gromacs use this to calculate the autocorrelation function
Screenshot_2020-12-08 Hydrogen bonds — GROMACS 2020-beta1 documentation
3. But which of my value 1. “Ac\sfin sys\v{}\z{}(t)” or 2. “Ac(t)” is this properly autocorrelation function?

When I have a function, how to compute hydrogen bond lifetime. I see that I can do this using this from instruction “with si(t)={0,1} for H-bond i at time t. The integral of C(τ) gives a rough estimate of the average H-bond lifetime τHB:”
Screenshot_2020-12-08 Hydrogen bonds — GROMACS 2020-beta1 documentation(1)
So first I should try obtain function form (but exponential will be e^x?) and then I definite integral from 0 to 5000 (my time is from 0 to 5000 ps)???

I read also something like that
“In general, these decays involve different simultaneous
processes, and hence the corresponding plots can be fitted to a
sum of exponentials”
Screenshot_2020-12-08 Sci-Hub Water Isotope Effect on the Phosphatidylcholine Bilayer Properties A Molecular Dynamics Simul...
where T i is the time constant and Ai the amplitude of the ith
individual decay process and N is the number of exponentials.
N is generally not known beforehand. Because of the nonorthogonality of the exponentials, unambiguous determination of the
proper values for Ti , Ai , and N is not easy."

  1. So how to calculate A and T?

Thanks in advance

Hi Jakub,

  1. The first plot (Ac_{fin,sys} is a tail-corrected autocorrelation function. Autocorrelation tails can be noisy and for longer time-scales move below and above the x-axis as shown in the upper image here

The first column shows an attempt to correct the autocorrelation function tail behaviour, the second column shows the raw data for the autocorrelation function. The third column shows the correlation between contacts and hydrogen bonds, telling you how many contacts lead to a hydrogen bond in your two analysed groups. The fourth plot is a simple numerical derivative of the first column.

  1. In your specific case, the negative values here are due to the tail correction, but generally autocorrelation functions exhibit the behaviour shown above, where for long time-scales curves decay in a sinoidal pattern, showing also negative values.

  2. The second column is the “true” autocorrelation function, the first one the corrected one Erik Marklund thought you should use to analyse time-scales. Fitting autocorrelation functions with noisy tails is a science in itself, for in-depth reading have a look here.

  3. To determine life-times, you would like to match an exponential function where you can read off the life-time from the exponent \tau from \exp^{-\frac{t}{\tau}} to the autocorrelation curve that you obtained.

One approach for stable fitting exponential decays with single exponentials is to match the areas of an ideal to the area under your autocorrelation function, cutting off noise tails or correcting for them. Some other ideas for calculating the life-time you can have a look here:

1 Like

Excuse me. I have a few questions about the third column. Does the “correlation” in the third column mean the “cross-correlation”? And how to evaluate the results? Does positive value in the correlation (or cross-correlation) mean more related between two analyzed group? How about the negative value? I wonder how the value was calculated.