Magnetic atoms serve as the ultimate building blocks at the atomic scale for nextgeneration quantum devices, with applications in high-density data storage, quantum computing, and sensing. The design of atomic quantum devices begins with surface adsorption, which enables single-atom addressability-the ability to initialize, control, and detect individual spins-while also allowing scalability to surface spin arrays. Recent advancements in scanning tunneling microscopy (STM) combined with electron paramagnetic resonance (EPR) spectroscopy have made it possible to control and detect individual atomic spins on surfaces. In the pursuit of scalability, STM techniques now enable the bottom-up, atom-by-atom construction and characterization of atomic spin lattices on surfaces. Currently available theoretical approaches for magnetic adsorbates lag behind these experimental breakthroughs. Density functional theory (DFT) with periodic boundary conditions (pbc) remains the workhorse for ab initio simulations of magnetic materials. However, DFT often fails to provide even a qualitatively correct description of key ground-state properties, such as adsorption energies and charge-transfer processes, which are crucial at adsorbate-surface interfaces. To compute magnetic properties, the Hubbard (U ) correction to DFT is widely used. However, choosing the appropriate U parameter can be challenging, often resulting in only qualitative agreement with experiment, or failing to predict correct state energy splitting and fundamental properties. While more accurate post-Hartree-Fock methods for periodic systems are available, they are mostly used for ground-state properties of dynamically-correlated materials and become computationally prohibitive as the number of atoms increases in the simulation cell. For local phenomena, such as molecular and atomic adsorbates on surfaces, quantum embedding theories provide a way to reduce computational effort without compromising accuracy. Among these methods, a groundbreaking development was the introduction of a density-based embedding scheme that combines multireference methods with periodic DFT for molecular adsorbates on surfaces. In this work, we adopt a more recent, projectionbased version of density embedding, which enforces orthogonality between fragment orbitals, eliminating the need for non-additive kinetic potentials or optimized embedding potentials that were required in previous formulations. Projection-based density embedding that combines equation-of-motion coupled-cluster singles and doubles (EOM-CCSD) theory with DFT has gained popularity for exploring both ground and excited states of molecular systems. Additionally, EOM-CCSD-in-DFT has been extended to open-shell species and was recently applied to investigate the spin states and magnetic properties of large transition-metal molecular magnets beyond the reach of standalone EOM-CCSD. Furthermore, some of us have applied embedded EOM-CCSD to elucidate dissociative electron attachment on metal surfaces, highlighting its effectiveness in addressing challenges in surface science. EOM-CCSD methods provide a robust and efficient framework for treating molecular magnets, describing their multiconfigurational wave functions within a single-reference approach. These methods have been extensively benchmarked across various transition-metal complexes, accurately predicting spin-state energy gaps, exchange couplings, g-tensors, spininversion energy barriers, and macroscopic magnetic properties such as magnetization and susceptibility. Additionally, the second-order approximate coupled-cluster singles and doubles (CC2) method has recently been applied to molecular magnets, accurately describing their spin states. However, the steep polynomial scaling of EOM-CCSD (N 6 ) and CC2 (N 5 ), where N is the system size, limits their applicability to small-to-medium-sized molecular systems. In this work, we employ the more cost-efficient embedded EOM-CCSD method to investigate the spin states and magnetic behavior of individual magnetic atoms on surfaces, going beyond traditional DFT approaches. As a case study, we focus on single Co atoms adsorbed on the MgO(001) surface. This system serves as an ideal platform for designing atomic-scale quantum devices. A combined experimental and DFT+U study has predicted that Co atoms on MgO(001) exhibit the highest magnetic anisotropy energy (or spin-inversion barrier) among 3d transition-metal adatoms, which is an essential feature for realizing efficient single-atom magnets. This behavior is attributed to the axial ligand field of Co at the oxygen adsorption site, which preserves the orbital angular momentum of gas-phase Co atoms (L = 3), enhancing spin-orbit coupling (SOC) and resulting in a record magnetic anisotropy energy of 58 meV (468 cm -1 ). In addition, the Co/MgO(001) system provides a simple yet insightful model for assessing the performance of EOM-CCSD-in-DFT compared to EOM-CCSD and CC2.