The bottom two layers were fixed to simulate the bulk and the top two layers were free to relax to resemble the surface. 15 Å of vacuum is added in the z direction (perpendicular to the surface) to avoid interaction between periodic images. The Monkhorst-Pack scheme was used for K-points of 4 × 4 × 1. The energy cut-off and the relative cut-off used were 550 and 50 Rydberg, respectively. The force convergence was taken to be 3 × 10 Bohr -1 Hartree. The exchange correlation functional of Perdew, Burke, and Ernzerhof (PBE) was used. All these parameters were chosen after running convergence and sensitivity tests. DFT calculations are all performed using CP2K code and further computational details are given elsewhere. CO was used as a starting intermediate as it one featured in the CO2 reduction pathway to higher order products. To compute energy differences of elementary proton coupled electron transfer (PCET) steps, the computational hydrogen electrode (CHE) model 37 was used. In this model we assume hydrogen gas is at equilibrium with proton and electron and the corresponding potential is 0 V vs. RHE, thus the energy of proton coupled with electron is estimated by half of the energy of hydrogen gas: