  The O vacancy (Ov) formation energy, $E_\textrm{Ov}$, is an important
property of a metal-oxide, governing its performance in applications such as
fuel cells or heterogeneous catalysis. These defects are routinely studied with
density functional theory (DFT). However, it is well-recognized that standard
DFT formulations (e.g. the generalized gradient approximation) are insufficient
for modeling the Ov, requiring higher levels of theory. The embedded cluster
method offers a promising approach to compute $E_\textrm{Ov}$ accurately,
giving access to all electronic structure methods. Central to this approach is
the construction of quantum(-mechanically treated) clusters placed within
suitable embedding environments. Unfortunately, current approaches to
constructing the quantum clusters either require large system sizes, preventing
application of high-level methods, or require significant manual input,
preventing investigations of multiple systems simultaneously. In this work, we
present a systematic and general quantum cluster design protocol that can
determine small converged quantum clusters for studying the Ov in metal-oxides
with accurate methods such as local coupled cluster with single, double and
perturbative triple excitations [CCSD(T)]. We apply this protocol to study the
Ov in the bulk and surface planes of rutile TiO2 and rocksalt MgO, producing
the first accurate and well-converged determinations of $E_\textrm{Ov}$ with
this method. These reference values are used to benchmark exchange-correlation
functionals in DFT and we find that all studied functionals underestimate
$E_\textrm{Ov}$, with the average error decreasing along the rungs of Jacob's
ladder. This protocol is automatable for high-throughput calculations and can
be generalized to study other point defects or adsorbates.
