  Quantum nuclear effects and anharmonicity impact a wide range of functional
materials and their properties. One of the most powerful techniques to model
these effects is the Stochastic Self-Consistent Harmonic Approximation (SSCHA).
Unfortunately, the SSCHA is extremely computationally expensive, prohibiting
its routine use. We propose a protocol that pairs machine learning interatomic
potentials, which can be tailored for the system at hand via active learning,
with the SSCHA. Our method leverages an upscaling procedure that allows for the
treatment of supercells of up to thousands of atoms with practically minimal
computational effort. The protocol is applied to PdCuH$_x$ ($x = 0-2$)
compounds, chosen because previous experimental studies have reported
superconducting critical temperatures, $T_\text{c}$s, as high as 17~K at
ambient pressures in an unknown hydrogenated PdCu phase. We identify a $P4/mmm$
PdCuH$_2$ structure, which is shown to be dynamically stable only upon the
inclusion of quantum fluctuations, as being a key contributor to the measured
superconductivity. For this system, our methodology is able to reduce the
computational expense for the SSCHA calculations by $\sim$96\%. The proposed
protocol opens the door towards the routine inclusion of quantum nuclear motion
and anharmonicity in materials discovery.
