We also applied this correction throughout the work. It is worth to mention practical limitations of N V E MD simulations due to the limited total number of collected frames N , determined by available computational and time resources. We calculate N t trajectories each with N s steps, so that the total number of frames is N = N t • N s . For the predefined constant N we may balance N t and N s but cannot increase both parameters. Smaller N s cause the lowering of the spectral resolution, because of the duration of trajectories t = N s • ∆t, where ∆t is the integration time step. On the other side, smaller N t will result in poorer phase-space sampling. To overcome this problem of N V E MD, we can try switching to N V T MD with a thermostat, in which the phase space sampling is done on-the-fly. In this case, with a similar number of computational operations it is possible to obtain fewer trajectories (smaller N t ) but with larger lengths (larger N s ) and optimal phase space sampling. However, the thermostat must not interfere too much with the dynamics of the system. Andersen thermostat is the most straightforward to implement. The WS method relies on equilibrium molecular geometry and cannot be used for anything except for sampling of initial conditions. SWS, on the other hand, has the same principle as MBS and thus can be used with the Andersen thermostat. Test calculations have been performed for difluoromethane CH 2 F 2 at the PBE/def2-SVP level of theory using a combination of SWS and Andersen thermostat (abbreviated here as SWS-AT+MD). We also did classical constant temperature MD simulations with the Andersen thermostat applying the MBS scheme (MBS-AT+MD). Figure shows the obtained vibrational spectra in comparison to harmonic and VPT2 fundamentals. Also, for demonstrational purposes, we calculated and plotted spectra from the WS+MD and SWS+MD simulations, which were used in computation of rotational constants described above. Classical MBS introduces relatively small vibrational amplitudes. Therefore, high-frequency vibrations from MBS-AT+MD are essentially harmonic. Taking into account nuclear quantum effects in WS+MD, SWS+MD, and SWS-AT+MD simulations leads to a better agreement of the MD-based spectra with the VPT2 calculations. Short lengths of MD trajectories in the N V E MD simulations limit the frequency resolution in the WS+MD and SWS+MD simulations. The much longer SWS-AT+MD simulations show results very similar to WS+MD and SWS+MD, which is seen in the similar shapes and positions of the peaks in the spectra. However, due to the much longer duration of trajectories The line around harmonic spectra is the Gaussian-convoluted spectra with each fundamental peak broadened with FWHM of 10 cm -1 . The lines above show the spectra obtained from the MD simulations.