We introduce a novel method for the simulation of the impact scattering in vibrational scanning transmission electron microscopy electron energy loss spectroscopy (STEM-EELS) simulations. The phonon-loss process is modeled by a combination of molecular dynamics and elastic multislice calculations within a modified frozen phonon approximation. The key idea is thereby to use a so-called δ-thermostat in the classical molecular dynamics simulation to generate frequency dependent configurations of the vibrating specimen's atomic structure. The method includes correlated motion of atoms and provides vibrational spectrum images at the cost comparable to standard frozen phonon calculations. We demonstrate good agreement of our method with simulations and experiments for a 15 nm flake of hexagonal boron-nitride (hBN).Excellent spatial resolution is the highlight of the transmission electron microscope [1-6], allowing atomically resolved elemental mapping and chemistry. Recent instrumental developments have improved the energy resolution of electron energy loss spectroscopy down to 4.2 meV [7,8]. This unique combination of high spatial and energy resolution opened doors to unprecedented experiments such as temperature measurement at the nanoscale [9, 10], identification and mapping of isotopically labeled molecules [11], position-and momentum-resolved mapping of phonon modes [12,13], mapping of bulk and surface modes of nanocubes [14], or investigations of the nature of polariton modes in van der Waals crystals [15].Inelastic electron scattering on atomic vibrations in a polar material consists of two major contributions, namely, impact scattering and dipolar scattering [16,17]. The latter stems from a long ranged interaction between the beam electron and an oscillating dipole moment generated by the atomic vibrations. It can be used to perform damage-free vibrational spectroscopy in an aloof beam geometry [18,19]. High-angle impact scattering on the other hand is quite localized scattering on the atomic potential and allows for high resolution imaging and spectroscopy [17,[20][21][22].From a theoretical point of view, there is a limited number of approaches one can follow to simulate vibrational spectroscopy and imaging. Dwyer presents a treatment of both single impact and single dipolar scattering including a detailed phonon dispersion from a density functional theory calculation in a multislice scheme [23]. Forbes et al. introduced a quantum excitation of phonons model, which allows for the inclusion of thermal diffuse scattering into multislice image simulations [24]. The thermal vibrations are modeled within the Einstein model for quick image calculations, neglecting correlations in atomic movements and dispersion effects. Such an approach is not able to simulate vibrational spectra but in principle one could use different models of atomic vibrations to generate the necessary configurations. In a more recent work, Forbes et al. described a method to model the vibrational sector of the electron energy loss spectrum [25...