Velocity autocorrelation function for molecular dynamics trajectory

Question

I have a 1 ns long molecular dynamics trajectory of 884 molecules of water, and I am trying to compute the velocity autocorrelation function using MDAnalysis to analyze the GROMACS output files.

I have written:

def vacf(ncorr):
    acf, norm = np.zeros(ncorr), np.zeros(ncorr) # initialize autocorrelation function and normalization arrays

    n = 0 # first frame
    for ts in u.trajectory: # loop through the trajectory
        vel = u.trajectory[n].velocities.flatten() # calculate velocities at t0
        maxn = min(u.trajectory.n_frames, n + ncorr) - n

        for i in range(maxn):
            veln = u.trajectory[n + i].velocities.flatten() # compute velocities at time t + t0
            acf[i] = acf[i] + np.dot(vel, veln.T) / (len(u.atoms))
            norm[i] = norm[i] + 1

        n += 1

    acf = acf / norm
    return acf

The velocity autocorrelation function results are usually something that start at 1 and slowly decay to 0, but this is not the case for the simulation I'm analyzing. The output I'm getting is:

Do you have any ideas on what I'm doing wrong? Thank you very much in advance, have a great day!

Answer 1

I think I managed to solve the problem. The code is now:

def vacf():
    acf, norm = np.zeros(int(ncorr)), np.zeros(int(ncorr))

    n = 0
    for atoms in u.atoms:
        vel = u.atoms[n].velocity
        maxn = min(u.trajectory.n_frames, n + ncorr) - n

        for i in range(maxn):
            veln = u.atoms[n + i].velocity
            acf[i] = acf[i] + np.dot(vel, veln.T) / (np.dot(vel, vel))
            norm[i] = norm[i] + 1

        n += 1

    acf = acf / norm
    
    return acf