A class of methods that addresses some of these challenges are selected CI methods, which truncate the size of the Fock space by selecting only those configurations that are relevant to describe the states of interest. In particular, the combined density functional theory and MRCI (DFT/MRCI) method has been shown to be able to achieve an accuracy comparable to CASPT2 with reduced computational cost, as well good performance for various spectroscopic applications. This method employs a selected CI algorithm for the generation of the reference space, which when combined with DFT-specific corrections to the Hamiltonian matrix, results in highly compact wave function expansions. In addition, our group has developed a perturbative approximate variant of the parent DFT/MRCI method, which is termed DFT/MRCI(2) and drastically reduces the computational cost while maintaining the accuracy of the original method. In this work, we will primarily employ DFT/MRCI (2) and expect the results and conclusions to be equally applicable to DFT/MRCI. While the generation of the reference space using selected CI is a core-feature and key strength of the DFT/MRCI(2) method, it also creates challenges in the construction of smooth PESs. In common with other selected CI methods, the general, non-full CI configuration space generated via the selection procedure renders the corresponding adiabatic energies non-invariant with respect to orbital rotations, complicating the implementation of analytical gradients. More importantly, the reference space and first order interaction space generation procedure, as implemented in DFT/MRCI (2), results in PESs that are locally non-smooth when the spaces generated are specific to each nuclear configurations. That is, it cannot be guaranteed that the DFT/MRCI(2) energies of a set of closely distanced structures constitute a reliably differentiable PES. However, the DFT/MRCI(2) PESs resemble continuous underlying surfaces that have been disrupted by deterministic noise that arises from the native selection procedure of the method, making it locally non-smooth. As a result, implementations of standard numerical gradients are not guaranteed to yield physical or meaningful quantities and will be strongly dependent on the step size. A natural way to address the aforementioned problem is by building a smooth surrogate representation of the PES, on which standard numerical or analytical gradient schemes can be employed. To this end, Gaussian process regression (GPR) stands out for its simplicity and versatility. First, GPR performs well with small amounts of training data and is less data intensive compared to neural network architectures, thus making it suitable for on-the-fly generation of PESs.