We explore the capabilities of the frequency-resolved frozen phonon multislice method introduced in Phys. Rev. Lett. 124, 025501 (2020) to model inelastic vibrational scattering in transmission electron microscopy. We review the method in detail and discuss advantages of using a so called hotspot thermostat instead of the $\delta$-thermostat used in our first report. We apply the method to simulate vibrational electron energy loss spectra of hexagonal boron nitride under plane wave illumination. Simulated spectroscopic information well represents the theoretical phonon bandstructure of the studied material, both in terms of energies as well as polarization vectors of individual phonon modes.