Ab initio molecular dynamics (AIMD) within the Born-Oppenheimer approximation is a promising alternative to overcome the limitations of the normal mode analysis, in particular for including anharmonicities along with temperature and solvent effects. However, this comes at the expense of increased computational cost in the case of systems with only a few hundreds of atoms, and the loss of nuclear quantum effects such as the vibrational zero-point energy or quantum tunneling when treating nuclei classically.