A longstanding challenge in computational heterogeneous catalysis is the inability of density functionals to accurately predict the adsorption energies and preferred adsorption sites of CO molecule on transition metal surfaces, a well-known "CO adsorption puzzle" . Experimental studies consistently show that at low coverages, where θ < 1/3 monolayer (ML), CO adsorbs on the top site of copper atoms . As the coverage increases, CO molecules begin to occupy twofold bridge and threefold hollow sites, indicating a coverage-dependent shift in adsorption behavior .