The Metropolis-Hastings Markov chain Monte Carlo algorithm was adapted for use in the Virtual Surface Site Relaxation-Monte Carlo (VSSR-MC) method . VSSR-MC simulations were conducted within the grand ensemble, where the chemical potentials of adsorbates were specified, allowing the total number of adsorbates to fluctuate throughout the MC run. In each iteration, a random adsorption site was selected for a potential state change. An intermediate oxygen chemical potential of μO = -2.5 eV was used to capture surface motifs relevant to both oxygen-poor and oxygen-rich conditions, and the structures were then re-weighted at other μO of interest. The initial sampling temperature was set to 12,000 K (kBT = 1). After each discrete sampling step, a continuous relaxation process was applied using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm within the ASE. To optimize computational efficiency and avoid divergence due to instability in the MLIP, the relaxation was limited to a maximum of 100 steps.