How exactly does one compute the vibrational density of states from the output of a molecular dynamics simulation?

Question

There are many things that can be done with the output of a molecular dynamics simulation, but one of the more powerful things is the ability to compute the vibrational density of states (or the infrared spectrum and raman spectrum). It is well known that the vibrational density of states can be related to the fourier transform of the velocity autocorrelation function (VACF). My question is, how exactly do I compute this autocorrelation function?

To be more concrete, suppose I have a large file of xyz-formatted velocities, as I would get from a molecular dynamics simulation (MD). Each velocity frame is separated in time by $k\Delta t$, where $k$ is the stride at which you save the velocities and $\Delta t$ is the actual simulation time step.

The VACF is defined as, $$ C(t)=A\langle\vec{v}(0)\cdot\vec{v}(t)\rangle $$ where $A$ is some normalization constant and $\vec{v}(t)$ is the velocity of the system at time $t$. The reason I am asking this question is not because there aren't resources explaining how to calculate this quantity, and why it is useful, but because there are way too many of them and they often contradict each other.

First of all, there is a fairly trivial way to calculate this which is to actually do the dot-product between all the velocity vectors for all values of time and just average them. This scales as $O(N^2)$ and is almost never used because there is an elegant way to calculate correlation functions using fourier transforms which is much more efficient, so the answer should describe that approach, though I don't think it's necessary to describe what a fourier transform actually is in answering this question.

Here are some questions I would like answered in detail, which aren't totally clear to me because of conflicting information I have seen.

• What value should be chosen for the normalization constant $A$?
• Is the final VACF meant to be the sum of the VACFs of the velocity components for each atom, or the mean of the VACFs? (I'm pretty sure it's the sum, but I've seen the mean written in a few places.)
• This is the question that is least clear to me. Where does the gap between frames, $k\Delta t$, enter the picture when calculating the VACF?

---

There is a very good outline of how to do this mathematically in this Physics SE question. What I am hoping for is something more like an algorithm or an answer written in pseudo-code.

Answer 1

Ok, I figured out the answer, so I will go ahead and answer my own question.

First, I will provide an overview in bullet points:

• [Optional] Generate the velocities at evenly spaced frames from the geometries by the equation $\frac{x(t)-x(t-\Delta T)}{n\Delta t}$ where $x(t)$ is the geometry at time $t$. $n$ is the lag between frames form the simulation that are saved and $\Delta t$ is the simulation time step.
• Read in the velocities generated in the previous step or which were saved from the simulation. Note that it is slightly better to use the ones from the simulation since they will tend to be slightly more accurate due to the fact that after the simulation one does not usually have access to all of the data which one does during the simulation
• Compute the average $\langle v(t)v(t+T)\rangle$ (more on this in a minute)
• Renormalize the above time correlation by its value at $t=0$, so that the max value is $1.0$
• Perform a fourier transform of the renormalized velocity autocorrelation function (VACF)
• Get the discrete frequencies used to build up the fourier transform for plotting purposes

Detailed Answer:

First, you need to get your hands on a bunch of velocities which are evenly spaced in time. Typically this time will not be the simulation time, but you will need to know what the time actually is.

Now, the original source of my confusion was due to the fact that nearly any paper you read which says they computed the VACF will use the following equation: $$\mathrm{VACF}=\langle v(0)v(t)\rangle$$ This is an incredibly terse notation, normally $\langle\rangle$ would imply an ensemble average which means you average some quantity for one or many atoms over time. The fact a time average and the true ensemble average are equal is due to the ergodic principle, which I won't describe here.

What you are computing in the VACF is really not an ensemble average though. You are averaging over all atoms and all time lags.

So, I will rewrite this equation as:

$$ \langle v(t)v(t+T)\rangle=\sum_{i=1}^{T}\sum_{j=1}^{N}v_j(t_i)\cdot v_j(t_i+n\Delta t) $$

Once again, $n\Delta t$ is the time-spacing between frames. $T$ is the total length of signal you have (i.e. number of points in the time-series). Notice that since each time windows will be getting shorter as $i$ increases, one usually performs the first sum up to $M/2$ as long as this is long enough that the correlation function has sufficiently decohered within that window. This ensures that one has the same statistical accuracy in each time window by keeping the length of each window the same. One can further increase the accuracy by taking a very long time series as one might get from MD and chopping it into independent segments and averaging the result of the above signal over those independent segments.

Finally, one usually renormalizes $\langle v(t)v(t+T)\rangle$ by $\langle v(0)v(0)\rangle$. That is, the first point in the TCF will equal $1$ and all other points should be bounded between $-1$ and $1$. This is proper to do since this is a correlation function, and correlations are usually normalized to lie between $-1$ and $1$.

Computational Details:

The pertinent equation I have written up above is quite easy to calculate numerically. It is just two nested for loops. However, this can be extremely slow to calculate as it scales quadratically, and one often has very long time-series and many atoms when doing MD. It turns out, for reasons I won't explain here, that the calculation of either autocorrelation functions or cross-correlation functions can be written in terms of a fourier transform and inverse fourier transform. In that case, rather than directly computing the product as above, one computes a fourier transform of the time series, takes the product of that series with itself, and inverse fourier transforms.

Getting the VDOS:

Going from a correlation function to something more physically meaningful is usually rather simple, as there are many physical observables which are directly related to some kind of TCF. So, to get the VDOS, which is what I happened to be asking about, one simple performs a fourier transform of the VACF. There is a last point which is that the fourier transform builds up a frequency-space representation of a time-domain signal from periodic basis functions (usually it's a complex basis formed of sines and cosines). So, to actually plot the VDOS, which is what you usually want, you need to get the frequencies of these basis functions and the corresponding intensities.

Practical Details:

If all you want are simple auto-correlations and cross-correlations, there is a small python package called tidynamics which can do this. It's also quite easy to implement these correlations calculations using either pyfftw or the numpy fft module. Note that to get the frequencies for the VDOS, you need the function np.fft.fftfreq.