  We devise a framework to compute accurate phonons in molecular crystals even
in case of strong quantum anharmonicity. Our approach is based on the
calculation of the static limit of the phononic Matsubara Green's function from
path integral molecular dynamics simulations. Our method enjoys a remarkably
low variance, which allows one to compute accurate phonon frequencies after a
few picoseconds of nuclear dynamics, and it is further stabilized by the use of
appropriate constrained displacement operators. We applied it to solid hydrogen
at high pressure. For phase III, our predicted infrared (IR) and Raman active
vibrons agree very well with experiments. We then characterize the crystalline
symmetry of phase IV by direct comparison with vibrational data and we
determine the character of its Raman and IR vibron peaks.
