Hi Tom and anyone who can help,
After some experimentation with the cpptraj autocorr and our own code we think that we have reconciled the differences between the two calculations.
I would be grateful if anyone can corroborate our observation and comment whether this is indeed what the code in amber is doing.
We noticed that short tracks of h-bonds that go rapidly to 0, fit very well to those that we calculate with our code. But tracks that decay over a long period of time, i.e., longer than about 300 ns, the auc calculated with amber had a larger slope with essentially the same shape as our auc. This led to the supposition that amber may be adding zeros at the end of the track to ensure that the limit of the track is going to zero at infinite time.
Indeed, when we added zeros at the end of the tracks and used our program the auc were identical.
Of course, that is not a proof that this is what is done in cpptraj.
So we would like to get a possible answer of how autocorr in cpptraj is doing the calculation.
As I wrote before we use the nocovar method because otherwise the auc becomes negative at long times.
Thanks,
Rami
On Aug 30, 2016, at 6:07 PM, Tom Kurtzman <simpleliquid.gmail.com<mailto:simpleliquid.gmail.com>> wrote:
Hi Roman, in response to your comment "In my opinion this is not possible
since the track is made of 0 and 1, so the autocorrelation cannot be < 0."
Often correlation functions (including autocorrelation functions) are
normalized so that their initial values are 1 and absence of correlation is
zero so instead of calculating
<A(t)A(t')> they actually calculate This is the same as calculating
<[A(t) -<A>][A(t')-<A>]> = <A(t)A(t')> - <A>^2. This allows for negative
numbers and ensures that the limit as t --> \infty is zero. Often the
function is then divided by its initial value to give an initial value
(upper bound) value of unity.
For the no-covariance calculation, I assume that the average values are not
subtracted and hence the limit and infinite time approaches <A>^2 which is
positive.
It is not unusual for autocorrelation functions to have negative values nor
to have long time non-zero tails (the velocity autocorrelation function is
a famous example of both of these) however, I'm not sure whether the long
time tail you are seeing is
1) real
2) fluctuations about zero due to sampling limitations (noise)
3) an artifact of an incorrectly calculated quantity.
Without seeing the data, I would initially suspect it's just noise.
Tom
On Mon, Aug 29, 2016 at 9:05 AM, Osman, Roman <roman.osman.mssm.edu<mailto:roman.osman.mssm.edu>> wrote:
Hello Amberites,
I was wondering whether anybody has used an autocorrelation of hydrogen
bonds.
Here is a cpptraj input that I ran:
# Read in trajectory
trajin mcmrjp_free_nw.dcd
hbond HB out mcmrjp_free_hb.dat :17,252,67,206,159,272,27,32,
228,245,262,267,165
,183,34,138,239,290,39,143,207,204,162 series uuseries mcmrjp_free_hb.dat
avgout
mcmrjp_free_avg.dat printatomnum
run
runanalysis lifetime HB[solutehb] out mcmrjp_free_lifet.dat
autocorr HB[solutehb] out mcmrjp_free_auc.dat lagmax 14500
run
The surprising thing is that I get in the aucorrelation function negative
numbers towards the end of the track.
In my opinion this is not possible since the track is made of 0 and 1, so
the autocorrelation cannot be < 0.
In fact if I run with
autocorr HB[solutehb] out mcmrjp_free_auc.dat lagmax 14500 nocovar
The autocorrelations become positive.
Am I doing something wrong?
Thanks for your help.
PS: I can provide the trajectory (large) and the output from the
autocorrelation.
Roman Osman
Professor of Structural and Chemical Biology
Mount Sinai School of Medicine
New York, NY 10029
(212) 659-8627
roman.osman.mssm.edu<mailto:roman.osman.mssm.edu><mailto:roman.osman.mssm.edu>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org<mailto:AMBER.ambermd.org>
https://urldefense.proofpoint.com/v2/url?u=http-3A__lists.ambermd.org_mailman_listinfo_amber&d=AwICAg&c=4R1YgkJNMyVWjMjneTwN5tJRn8m8VqTSNCjYLg1wNX4&r=80h8gfWnXRd-mXPSrw98JBRjhncHe1Sopk2U91C1sZg&m=-hoofXSvixmoA1iOwfXi41mSpf3DaEZs1kdxf6s0p-Q&s=YXoZ0b2W1iVWNm4t_-F4rjkMCoaEH-OAP15eidtQO2s&e=
--
************************************************
Tom Kurtzman, Ph.D.
Assistant Professor
Department of Chemistry
Lehman College, CUNY
250 Bedford Park Blvd. West
Bronx, New York 10468
718-960-8832
https://urldefense.proofpoint.com/v2/url?u=http-3A__www.lehman.edu_faculty_tkurtzman_&d=AwICAg&c=4R1YgkJNMyVWjMjneTwN5tJRn8m8VqTSNCjYLg1wNX4&r=80h8gfWnXRd-mXPSrw98JBRjhncHe1Sopk2U91C1sZg&m=-hoofXSvixmoA1iOwfXi41mSpf3DaEZs1kdxf6s0p-Q&s=o3uJDn8But0uJdUcDMpoKaCqj2NRM6dxD3CbDdwWrHM&e=
<https://urldefense.proofpoint.com/v2/url?u=http-3A__www.lehman.edu_faculty_tkurtzman_index.html&d=AwICAg&c=4R1YgkJNMyVWjMjneTwN5tJRn8m8VqTSNCjYLg1wNX4&r=80h8gfWnXRd-mXPSrw98JBRjhncHe1Sopk2U91C1sZg&m=-hoofXSvixmoA1iOwfXi41mSpf3DaEZs1kdxf6s0p-Q&s=G9qKDa1ESYL7tKFnPsa8zDrGteWJqxHWTUzqap0HM78&e= >
************************************************
_______________________________________________
AMBER mailing list
AMBER.ambermd.org<mailto:AMBER.ambermd.org>
https://urldefense.proofpoint.com/v2/url?u=http-3A__lists.ambermd.org_mailman_listinfo_amber&d=AwICAg&c=4R1YgkJNMyVWjMjneTwN5tJRn8m8VqTSNCjYLg1wNX4&r=80h8gfWnXRd-mXPSrw98JBRjhncHe1Sopk2U91C1sZg&m=-hoofXSvixmoA1iOwfXi41mSpf3DaEZs1kdxf6s0p-Q&s=YXoZ0b2W1iVWNm4t_-F4rjkMCoaEH-OAP15eidtQO2s&e=
Roman Osman
Professor of Structural and Chemical Biology
Mount Sinai School of Medicine
New York, NY 10029
(212) 659-8627
roman.osman.mssm.edu<mailto:roman.osman.mssm.edu>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Sep 01 2016 - 09:30:02 PDT
