I have not looked at this code in a while but I assume that correlation functions in cpptraj are done through the convolution theorem using fast fourier transforms. To do this efficiently, you need to pad the end of your time array with zeroes to get to a power of 2. This seems to agree what you are seeing.
adrian
Sent from my iPad
> On Sep 1, 2016, at 10:07 AM, Osman, Roman <roman.osman.mssm.edu> wrote:
>
> 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
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Sep 01 2016 - 09:30:03 PDT
