Dynamical analysis

One strength of molecular dynamics with respect to other methods such as Monte Carlo is the fact that the time evolution of the system is followed and--therefore--information on the dynamics of the system is fully present. Experimental data (for instance, phonon dispersion curves in the case of crystals) are usually available in the wavevector-frequency space. Information obtained by MD in real space and real time will therefore be Fourier-transformed with respect to both space and time for useful comparisons with experimental results. To perform dynamical analysis, one typically instructs the simulation program to dump periodically on a mass storage medium the positions and possibly the velocities of all the particles. These can easily produce very large files: the calculation of dynamical quantities requires plenty of room for data. These data are then analyzed later by other programs to extract the quantities of interest. In some cases, dynamical analysis are directly performed by the MD program while it is running, thus saving disk space.

One typical quantity is the Van Hove correlation function $G({\bf r},t)$, representing the probability to find an atom at position ${\bf r}$ at time t if there was an atom at position ${\bf r}=0$ at time t=0. This quantity is connected by a double Fourier transform to the dynamic structure factor $S({\bf k},\omega)$reported by experiments:  

$$S({\bf k},\omega) = \int\int d{\bf r} \, dt \, e^{i ({\bf k}\cdot{\bf r} - \omega t )} G({\bf r},t) .$$ (21)

As for static quantities, the most convenient way to proceed in terms of computational effort is that of using the space Fourier transform of the density as a building block:  

$$\rho ({\bf k},t) = \sum_i {\rm e}^{i {\bf k}\cdot{\bf r}_i(t)}$$ (22)

This is a single-particle quantity which is very fast to compute. From here, one constructs a time-dependent density-density correlation function (also called intermediate scattering function) by performing an average over a molecular dynamics trajectory:

$$F({\bf k},t) = \frac{1}{N} \left\langle \rho ({\bf k},t+t_\circ) \rho (-{\bf k},t_\circ) \right\rangle .$$ (23)

The average is made for all the available choices for $t_\circ$.The maximum t (connected with the frequency resolution) is clearly limited by the length of the run, while the periodicity with which data are collected--typically a few time steps--controls the maximum frequency in the final spectra. The dynamic structure factor is finally obtained by a time Fourier transform:  

$$S({\bf k},\omega) = \int dt\, {\rm e}^{-i\omega t} F({\bf k},t)$$ (24)

Again, the allowed ${\bf k}$ are quantized as a consequence of periodic boundary conditions. $S({\bf k},\omega)$ will exhibit peaks in the $({\bf k},\omega)$plane, corresponding to the dispersion of propagating modes. The broadening of these peaks is related with anharmonicity. By performing calculations at different temperatures, one can for instance study the effect of temperature on the phonon spectra, fully including anharmonic effects: in this sense, the approach is clearly superior to standard lattice dynamics techniques. The main inconvenience of MD is that it is not so straightforward to extract eigenvectors.

When dealing with vibrational modes and different polarizations, it is often convenient to follow the approach outlined above, but using currents  

$${\bf j} ({\bf k},t) = \sum_i {\bf v}_i(t) {\rm e}^{i {\bf k}\cdot{\bf r}_i(t)}$$ (25)

in place of the density $\rho ({\bf k},t)$.