Voids in face-centered cubic (fcc) metals are commonly assumed to form via the aggregation of vacancies; however, the mechanisms of vacancy clustering and diffusion are not fully understood. In this study, we use computational modeling to provide a detailed insight into the structures and formation energies of primary vacancy clusters, mechanisms and barriers for their migration in bulk copper, and how these properties are affected at simple grain boundaries. The calculations were carried out using embedded atom method (EAM) potentials and density functional theory (DFT) and employed the site-occupation disorder code (SOD), the activation relaxation technique nouveau (ARTn) and the knowledge led master code (KLMC). We investigate stable structures and migration paths and barriers for clusters of up to six vacancies. The migration of vacancy clusters occurs via hops of individual constituent vacancies with di-vacancies having a significantly smaller migration barrier than mono-vacancies and other clusters. This barrier is further reduced when di-vacancies interact with grain boundaries. This interaction leads to the formation of self-interstitial atoms and introduces significant changes into the boundary structure. Tetra-, penta-, and hexa-vacancy clusters exhibit increasingly complex migration paths and higher barriers than smaller clusters. Finally, a direct comparison with the DFT results shows that EAM can accurately describe the vacancy-induced relaxation effects in the Cu bulk and in grain boundaries. Significant discrepancies between the two methods were found in structures with a higher number of low-coordinated atoms, such as penta-vacancies and di-vacancy absortion by grain boundary. These results will be useful for modeling the mechanisms of diffusion of complex defect structures and provide further insights into the structural evolution of metal films under thermal and mechanical stress.