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$):
$$\mathbf{v}_n(\omega)=\int_{-\infty}^{\infty}\mathbf{v}_n(t)e^{i\omega{}t}dt$$
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.
