For the DFT/MRCI(2) MECI optimizations, we consider a distribution of the parameters g, h, s x , s y , P, B and r IJ computing a set of 50 surrogates, each constructed with different LHS centered at the same geometry. The reason for this is that we found that using atom-centered Cartesian LHS to generate the training data can result in deviations between different learned surrogate surfaces. This can then manifest itself in slightly dif-ferent optimised MECI geometries and branching space topographies. This dependence on the training data is expected due to GPR being a non-parametric method. In contrast to parametric fitting, GPR does not assume a functional form, but is based solely on the correlation between each training point, given by the choice of kernel, thus its sensitivity to the training set. We have found this to be a more general problem, and believe LHS generated nuclear configurations are not ideally suited to general, on-the-fly construction of PESs when the kernel hyperparameters are obtained from MLE, in terms of both consistency and scalability. In addition to the curse of dimensionality from which LHS suffers, it has been shown that GPR hyperparameters optimized with MLE from data generated with LHS can perform poorly due to irregular pairwise distance distribution. The question of optimal sampling and fitting, however, is not addressed in this work and we stick to atom-centered Cartesian LHS. To confirm the observed sensitivity does not instead originate from the underlying DFT/MRCI(2) calculations, however, the same procedure of constructing multiple different surrogates was employed using ab initio MRCI energies, and the same sensitivity of the parameters on the training data was found. This will be discussed in more detail subsequently.