Computational detailsCCSD(T) calculations of CO, CO 2 , and H 2 adsorption energies at Ni Mg defects at the MgO(001) surface, as well as formation energies of bulk Ni Mg , are performed using ORCA electronic-structure package, which is based on Gaussian-type orbitals [1]. For bulk Ni Mg defects, NiMg 12 O 6 (PP) 50 cluster is used (Fig. 1c in main ). 14 electrons are removed to ensure bulklike ionization states of all the atoms. The positions of the atoms, ECPs, and point charges in the bulk cluster models are generated by cutting out a cube centered on a Mg atom from the bulk MgO structure with 4.211 Å) lattice constant, with no further relaxation. The lattice constant was obtained from standard HSE06 calculations (see below for details). For surface Ni Mg defects, two models were used: terrace-site and monolayer-step ( Fig. 1a-b in main text), both have the stoichiometry NiMg 8 O 9 (PP) 49 . The structure of the terrace-site cluster was generated as follows. First, a periodic slab model of MgO (001) (2x2 surface unit cell and 4 atomic layers) with a Ni Mg defect and possibly an adsorbed molecule is relaxed at HSE06 level with tight basis set and 4x4x4 k-points. Second, the cluster around the Ni Mg defect was cut out of the slab and immersed into an array of point charges and pseudopotentials. The geometry of embedding corresponded to ideal bulk-like structure. Since the relaxation of (sub)surface atomic layers is pronounced very slightly even with Ni, such embedding geometry matches cluster structure well. In the case of monolayer-step cluster, similar strategy was applied. The 4-atomic layer MgO(001) slab (2x3 surface unit cell) with a periodically translated along x-axe (2x2) "monolayerisland" and with/without adsorbate molecule was fully relaxed at HSE06 level. The cluster with Ni Mg defect was cut from it and embedded like in the other cases. For all cluster models, the target properties were converged with respect to the number of shells of pseudopotentials and point charges at PBE level. Collinear spin-unrestricted calculations are performed, unless otherwise specified, with spinmoment projection fixed to S z = 1, which is found to be the ground spin-state of the neutral Ni Mg defect for all α. The spin-states of charged defects are discussed in the main text.In the CCSD(T) calculations, cc-pVXZ [3] basis sets are employed, with X = D for Mg and O, and X = T for all atoms in adsorbed species and surface atom underneath (Ni or O). For adsorption and defect formation energies, counterpoise basis-set superposition error (BSSE) correction is applied. The excitations from 1s-3p orbitals of Ni, 1s-2p of Mg, 1s of O and C were excluded from the coupledcluster calculations (frozen-core approximation). These settings were the best we could afford with our embedded cluster models. However, we estimate the complete basis-set (CBS) limit for CCSD(T) by combining the extrapolation of MP2 energies to CBS limit with the so-called "focal-point" method [4] which is based on the observation that the difference in correl...