  The study of molecular impurities in $para$-hydrogen ($p$H$\rm_2$) clusters
is key to push forward our understanding of intra- and intermolecular
interactions including their impact on the superfluid response of this bosonic
quantum solvent. This includes tagging with one or very few $p$H$\rm_2$, the
microsolvation regime, and matrix isolation. However, the fundamental coupling
between the bosonic $p$H$\rm_2$ environment and the (ro-)vibrational motion of
molecular impurities remains poorly understood. Quantum simulations can in
provide the necessary atomistic insight, but very accurate descriptions of the
involved interactions are required. Here, we present a data-driven approach for
the generation of $impurity\cdots p$H$\rm_2$ interaction potentials based on
machine learning techniques which retain the full flexibility of the impurity.
We employ the well-established adiabatic hindered rotor (AHR) averaging
technique to include the impact of the nuclear spin statistics on the
symmetry-allowed rotational quantum numbers of $p$H$\rm_2$. Embedding this
averaging procedure within the high-dimensional neural network potential (NNP)
framework enables the generation of highly-accurate AHR-averaged NNPs at
coupled cluster accuracy, namely CCSD(T$^*$)-F12a/aVTZcp in an automated
manner. We apply this methodology to the water and protonated water molecules,
as representative cases for quasi-rigid and highly-flexible molecules
respectively, and obtain AHR-averaged NNPs that reliably describe the H$\rm
_2$O$\cdots p$H$\rm_2$ and H$\rm _3$O$^+\cdots p$H$\rm_2$ interactions. Using
path integral simulations we show for the hydronium cation that umbrella-like
tunneling inversion has a strong impact on the first and second $p$H$\rm_2$
microsolvation shells. The data-driven nature of our protocol opens the door to
the study of bosonic $p$H$\rm_2$ quantum solvation for a wide range of embedded
impurities.
