We set out to iteratively construct a versatile MLFF, which can subsequently be used in combination with an evolutionary search, specifically the covariance matrix adaptation evolution strategy (CMA-ES), to explore different unit cells of the Ti-rich SrTiO 3 (110) surface. Due to the complex nature of the various surface reconstructions, the MLFF must be capable of resolving a wide range of local environments. Additionally, given the nature of the ES, which produces more diverse and possibly unphysical structures than, e.g., molecular dynamics, intermediate configurations are likely to exhibit unusual properties, such as unphysical bond lengths. For that reason, the MLFF needs to be robust and the underlying training data diverse, making the construction of a suitable training set containing ab-initio energies and forces far from trivial. Furthermore, we aim to explore (2 × 5) surface reconstructions with up to 450 atoms per unit cell where the computational costs associated with DFT evaluation is prohibitive for constructing diverse training databases. The MLFF must therefore also be able to generalize from smaller training structures to the larger unit cells explored. We started from the basis of our previous results on a Ti-poor reconstruction of SrTiO 3 (110), namely the (4 × 1) (see Fig. ). We constructed a database of 495 structures, re-evaluated using VASP with the r 2 SCAN functional and trained a ten-member committee based on the descriptor-based NeuralIL architecture. In committees the uncertainty is approximated by training a set of models that vary by initialization seed, hyperparameters, architecture, or training data, and monitoring their disagreement on a prediction to obtain the model variance. The majority of computational cost incurred when training a descriptorbased model can be attributed to the calculation of the descriptors and the associated vector-Jacobian product operator. NeuralIL enables a particularly efficient committee implementation by reusing these elements for all members, so that committees needed for uncertainty estimation can be trained with a negligible performance penalty. We then generated CMA-ES trajectories for the (2 × 2), (2 × 3)a and (2 × 3)b recon-structions. The CMA-ES samples a population of λ individuals x g k , k = 1, .., λ, for every generation g from the multivariate normal distribution