Our work covers the following three main aspects: 1. We show how to include expensive HiFi calculations within a MuFi approach to get a practically usable and accurate PESs for N ×1DHR treatments of larger molecular systems. The approach is put to the test in this study for systems of broad scientific and industrial relevance. 2. We discuss the limitations of the approach: since it is data-driven, the model can never be better than the underlying data. Still, the MuFi model provides a useful PES and uncertainty measures. 3. We explore how to use linear GPR MuFi models with Bayesian methods to increase the effectiveness of sampling. This allows to optimally predict locations of further calculations and thereby saves computational time. The prerequisite for this is a problem-oriented acquisition function whose maximum will determine the next calculation. We adapt the Bayesian method because there is no optimisation problem per se, and develop a new acquisition function type that takes several objectives into account based on thermodynamic principles.