Here, we present a general framework for computing the infrared anharmonic vibrational spectra of polyatomic molecules using Born–Oppenheimer molecular dynamics (BOMD) with PyRAMD software. To account for nuclear quantum effects, we suggest using a simplified Wigner sampling (SWS) approach simultaneously coupled with Andersen and Berendsen thermostats. We propose a new criterion for selecting the parameter of the SWS based on the molecules’ harmonic vibrational frequencies and usage of the large-time-step blue shift correction, allowing for a decrease in computational expenses. For the Fourier transform of the dipole moment autocorrelation function, we propose using the regularized least-squares analysis, which allows us to obtain higher-frequency resolution than with the direct application of fast Fourier transform. Finally, we suggest the usage of the pre-parameterized scaling factors for the IR spectra from BOMD, also providing the scaling factors for the spectra at the BLYP-D3(BJ)/6-31G, PBE-D3(BJ)/6-31G, and PBEh-3c levels of theory.