All calculations were performed using density functional theory (DFT) implemented in the plane-wave based Vienna ab initio simulation package (VASP). Electron-core interactions were described by the projector augmented wave (PAW) method. The generalised gradient approximation (GGA) according to Perdew, Burke, and Ernzerhof (PBE) was employed as exchange-correlation functional. In order to describe the magnetic properties of Cu 2 O and CuO, spin polarization was included for all cases, and the Hubbard parameter, U, was applied to all Cu 3d electrons according to Dudarev et al. Using this PBE+U method, the lattice parameter, energy of formation, magnetic moment, and band gap of Cu 2 O and CuO were investigated for various values of the Hubbard U parameter, as shown in the supporting information (Table ), and compared to experimental data. Best matching results were obtained for a U value of 8 eV, and therefore this value was used for all subsequent calculations. The same value was obtained by Zivkovic et al., who also performed a systematic study to determine the U value for copper oxides. For the Cu 2 O and CuO bulk phase, a planewave cutoff energy of 550 eV was applied using Gaussian smearing with a width of 0.05 eV. The Brillouin zone was sampled using a 12×12×12 Monkhorst-Pack 34 k -point grid for Cu 2 O and a 10 × 14 × 10 grid for CuO. The structure relaxations were assumed to be converged when the energy differences of the electronic self-consistent field (SCF) were smaller than 10 -5 eV and all forces were at least smaller than 1.0 meV Å-1 . Based on the optimized bulk structures, symmetric slabs with two identical surfaces were modeled for all possible terminations of the lowindex (100), (110), and (111) surfaces of Cu 2 O and CuO. For all surfaces, the surface energy was calculated for slabs with different thicknesses and was assumed to be converged when the difference in surface energy was smaller than 1.0 meV Å-2 . Based on these results, we chose a slab thickness between 17 Å and 21 Å, depending on the respective surface. In order to avoid interactions between the periodic slabs in z direction, a vacuum region of 20 Å was added in all cases. For the surface calculations, the plane-wave cutoff energy was set to 500 eV, and Gaussian smearing was applied with a width of 0.05 eV. The atomic position of every atom was fully relaxed until the Hellmann-Feynman forces were smaller than 1.0 meV Å-1 and the energy difference of consecutive SCF cycles was smaller than 10 -5 eV. The Brillouin zone was sampled using Γ-centered Monkhorst-Pack 34 k -point grids of 9 × 9 × 1, 9 × 7 × 1, 7 × 7 × 1, 11×7×1, 7×7×1 and 7×7×1 for Cu 2 O(100), Cu 2 O(110), Cu 2 O(111), CuO(100), CuO(110) and CuO(111), respectively. To achieve high precision density of states (DOS) calculations, the convergence criterion for the SCF cycle was decreased to 10 -8 eV and the density of the k -point mesh was increased significantly (at least doubled). The VASP-KIT package was used to post-process the DOS calculations. The charge on atoms was calculated using the Bader charge analysis and crystal structure drawings were visualized using VESTA.