How to calculate the velocity autocorrelation function in a molecular dynamics simulation?

Question

I am running a 3D molecular dynamics solution of a Lennard-Jones fluid with $N$ particles in a box, and I want to plot the velocity autocorrelation function. From this resource, I see that $$\psi(t) = \left\langle \sum_{i=1} ^N v_i(t) \cdot v_i(0) \right\rangle$$

is the non-normalized VACF. My problem is that I don't understand how to evaluate such a thing for my simulation.

The angle bracket convention from statistical mechanics is a convention used for time-averages. So, I know that $$\langle A \rangle = \lim_{\tau \rightarrow \infty }\int_0^{\tau} A(t) dt$$ or $$\langle A \rangle = \lim_{N\rightarrow \infty}\frac{1}{N}\sum_{i=1}^N A_i$$

So for my simulation which I run once, how do I find the average value of $$\sum_{i=1} ^N v_i(t) \cdot v_i(0) \quad (??)$$

This is how I run my simulation:

1. Initialize positions and velocities of particles
2. Evaluate forces on particles (using $F=-\nabla U$)
3. Perform a time-step using a symplectic integrator with some time-step $\Delta t$
   • find positions and velocities of particles in the new configuration
4. Repeat steps 2-3 until you reach your end time

I could run this simulation multiple times, and then take the velocity averages, but if my initial conditions and parameters are the same, I should get the same results... Any advice would be appreciated.

Answer 1

Since it is expected that your velocities will be randomised in the beginning (the VACF is an ensemble average, after all), you shouldn't get the same results when repeating the experiment. However, there is a much easier and more efficient way to evaluate the VACF accurately within a single simulation: you evaluate the autocorrelation function every $\Delta t$ steps. This can be done, because the VACF should be unchanged with respect to the absolute time. For example, you calculate $\pmb{v}(t) \cdot \pmb{v}(0)$ and then $\pmb{v}(t + \Delta t) \cdot \pmb{v}(\Delta t)$ up to your final timestep. If you take the average of these, you should arrive at the converged ensemble average fairly quickly. That is to say:

$$\psi(t) \approx \frac{1}{N_{steps}-\frac{t}{\Delta t} + 1}\sum_{i=0}^{N_{steps}-t/\Delta t} \pmb{v}(t+i\Delta t) \cdot \pmb{v}(i\Delta t) $$

Finally, I would recommend to check out the fairly recent Matter Modeling Stack Exchange, where these types of questions are more likely to be answered.