Autocorrelation Function

[horlust]

Hello everyone!

I have calculated the autocorrelation function for the hydrogen bond lifetime in a system containing 884 molecules of water according to Gower's paper with MDAnalysis. The autocorrelation function is defined as:

C(\tau) = \langle \frac{h(t_0)h(t + t_0)}{h(t_0)^2} \rangle

However, to study the kinetics of the hydrogen bond I wanted to slightly modify the function according to Luzar's papers:

c(t) = \langle h(0) h(t) \rangle / \langle h \rangle

So I went to see the autocorrelation function in MDAnalysis' documents and I have:

# calculate autocorrelation
    for t in range(0, len(list_of_sets), window_step):
        Nt = len(list_of_sets[t])

        if Nt == 0:
            continue

        for tau in tau_timeseries:
            if t + tau >= len(list_of_sets):
                break

            # continuous: IDs that survive from t to t + tau and at every frame in between
            Ntau = len(set.intersection(*list_of_sets[t:t + tau + 1]))
            timeseries_data[tau - 1].append(Ntau / float(Nt))

    timeseries = [np.mean(x) for x in timeseries_data]

    # at time 0 the value has to be one
    tau_timeseries.insert(0, 0)
    timeseries.insert(0, 1)

Perhaps the person who wrote this function could clarify me this. From what I understand, list_of_sets is a list of the hydrogen bonds in each frame, so Nt is the number of hydrogen bonds in each frame. Ntau is the number of hydrogen bonds that still exist in a time t0 + t, then the average value is taken. It seems logical to me, but I was hoping to somehow find the direct calculation of the h(t) values, being it 1 if the hydrogen bond still exists and 0 otherwise. So I did something similar to what I thought to analyze the output:

def correlation_test(list_of_sets, tau_max=40, window_step=1):
    frames = len(list_of_sets)
    tau_timeseries = list(range(tau_max))  # start at 0
    timeseries_data = [[] for _ in range(tau_max)]

    # calculate correlation function
    for t in range(0, frames, window_step):
        for pair_set in list_of_sets[t]:
            h_t0 = 1 if pair_set else 0  # h(0)

            for tau in tau_timeseries:
                if t + tau >= frames:
                    break
                
                # calculate correlation function for each pair of frames in the window
                h_t = 1 if pair_set in list_of_sets[t + tau] else 0  # 1 if common pair, else 0
                timeseries_data[tau].append(h_t / h_t0)
    
    # calculate the average correlation function over all windows for each tau
    average_timeseries = [np.mean(x) for x in timeseries_data]
    
    # at time 0, the value has to be one
    tau_timeseries.insert(0, 0)
    average_timeseries.insert(0, 1)

    return tau_timeseries, average_timeseries, timeseries_data

The output is not similar to the one calculated by MDAnalysis (I've attached the image).

So I was hoping to get a better understanding of the autocorrelation function if possible! I was also hoping to calculate the autocorrelation function as described by the works of Luzar, maybe this could be implemented on MDAnalysis someday?

Thank you very much!
Best,
Maria