# Velocity autocorrelation function for molecular dynamics trajectory

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!

• Why code it by hand when you can just use gmx velacc (which is likely much faster than your Python implementation, less likely to contain bugs, and has more features)? Commented Feb 19 at 22:06
• Short answer: I didn't know it existed, thanks. Long answer: I like coding my own functions just to see if it works :) Commented Feb 29 at 8:04

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