I'm using molecular dynamics and I autocorrelate the velocities and Fourier transform them to obtain the phonon density of states (DOS). I have many doubts about this:

  1. The definition of DOS is: number of states per interval of energy at each energy level that ARE AVAILABLE to be occupied. If this is ok, then what I obtain autocorrelating velocities shouldn't be the DOS. I think that what I obtain is the population of the states, because with the velocities that the system have you should obtain the phonons that the system have, not the phonons that the system could have. I think it is something like the product of the probability distribution function and the DOS, but not just the DOS. Am I wrong? And if what I think is correct, the probability distribution function is the Bose-Einstein distribution (because the phonons are quantum particles) or the Maxwell-Boltzmann distribution (because of the classic behaviour of the atoms in my simulations)?

  2. In all places that I read about this they say that the system is in equilibrium. If I do a non-equilibrium molecular dynamic (NEMD) simulation, can't I correlate the velocities and obtain information about the phonons? Can I take a small part of the system, assume that this small part is in local equilibrium, correlate the velocities of the atoms in this part and obtain information about the phonons? What would I obtain?


1 Answer 1


I recently stumbled on a more or less similar question. The Fourier transform of the velocity autocorrelation does not give you the phonon DOS but rather the phonon population of your system. In other words, it gives you normal modes.

The actual calculation of the DOS can be found, for example, in the thesis of Dr. Hugo Ruiz, gracefully linked by Prof. G. Naumis in this discussion: https://www.researchgate.net/post/How_do_I_calculate_Phonon_Density_of_states_from_VACF

As Prof. Naumis says, the actual calculation is rare in the literature and it could clarify some questions to develop it here.

Unfortunately, it is in Spanish but I can give here a quick explanation of his calculations.

A. Computation of the phonon DOS

Given the Fourier transform of the velocities from the space of time ($t$) to the space of frequencies ($\omega$):


where the subscript $n$ corresponds to the $n^{th}$ atom and $i$ is the imaginary unit. From here one can get the spectrum of potential kinetic energy:

$$|\mathbf{v}_n(\omega)|^2 = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \mathbf{v}_n(t') \mathbf{v}_n(t) e^{i\omega{}(t-t')}dtdt'$$

In a stationary oscillatory state, the coordinate $r_n(t)$ can be rewritten as a function of the normal modes of oscillation:

$$ r_{nj}(t) = \sum_s Q_{snj} e^{-i\omega_{sj} t} $$

with the $\omega_s$ the normal frequencies of oscillation and j the direction of the 3D space, $Q_{snj}$ being the mean coordinate of the particle. Hence one can express velocities as the derivative regarding time:

$$ v_{nj}(t) = \sum_s Q_{snj}(-i\omega_{sj}) e^{-i\omega_{sj} t} $$

Using this expression in the second integral, one can deduce:

$$ \sum_{i=1}^N |\mathbf{v}_n(\omega_s)|^2 = \sum_s\sum_{(n*j)=1}^{3N} \int_{-\infty}^{\infty} |Q_{snj}|^2\omega_s^2 e^{i(\omega+\omega_s)t''}dt'' $$

where $t''=t-t'$.

At thermal equilibrium, given the equipartition of the energy, one has directly $\sum_j |Q_{snj}|^2\omega_{sj}^2=3k_BT$. As the normal modes are supposed to be dominant, the integral of the complex exponential can be reduced to a Dirac distribution such as $\rho(\omega) = \sum_s \delta(\omega+\omega_s)$ is the phonon DOS.

Hence: $$ \sum_{i=1}^N |\mathbf{v}_n(\omega_s)|^2 = 3Nk_BT\sum_s\delta(\omega+\omega_s) $$ $$ \downarrow $$ $$ \rho(\omega) = \frac{\sum_{i=1}^N |\mathbf{v}_n(\omega_s)|^2}{3Nk_BT} $$

Taking $t'=0$ for the origin of the time, this final expression reduces to $$ \rho(\omega) = \frac{\int_{-\infty}^{\infty}\sum_{i=1}^N <\mathbf{v}_i(t)\cdot\mathbf{v}_i(0)>}{3Nk_BT} $$ linking it to autocorrelation. If you just Fourier transform the autocorrelation of the velocities, you get the normal modes. You have to have them (i.e integrate) to get the actual DOS. More information on the assumptions is in the thesis.

B. Concerning phonons

Phonons follow a Bose-Einstein distribution as they can be created and annihilated by energy fluctuations just like bosons particles (photons for example) independently of a classical/quantum behavior of the system. As such they follow a Bose-Einstein distribution with a chemical potential equal to zero:

$$ <n_i> = \frac{1}{e^{\frac{\hbar\omega_i}{k_BT}-1}} $$

Furthermore concerning NEMD it depends on the kind of constraint you impose on your system. As you can see in the derivation of the DOS, the equipartition of the energy and thermal equilibrium is an important assumption. The main problem in NEMD will be the fluxes, not only of matter (which can happen even in solids and I don't know what kind of system you are studying) but also of heat which can end up causing some flux of momentum. The problem with fluxes is that they break the isotropy of your system and you can end up with local non-equipartition of the kinetic energy, breaking one of the assumptions of the derivation. To assume local equilibrium you have to be sure that the input and output of energy and/or matter compensate in every direction for a long enough time and on a large enough region to do your statistics. This is far to be obvious even (or especially) in stationary out-of-equilibrium states.

I think a better way to determine the normal modes of higher energy is to simulate your system at a higher temperature in the canonical ensemble with a good thermostat (Nosé-Hoover Chain if you can).

I hope this helps. Any comment is welcome.


Not the answer you're looking for? Browse other questions tagged or ask your own question.