Among the different phases of TiO2, the rutile form has been extensively studied, especially under reductive conditions, where surface reconstructions occur. Both experimental studies and theoretical simulations, such as density functional theory (DFT) calculations, have identified many stable surface terminations . However, capturing the full range of surface configurations remains computationally challenging due to the vast number of possible structures that need to be sampled. The primary challenge lies in identifying nontrivial, complex surface reconstructions that significantly influence the material's catalytic behavior.