Cobalt nanoparticles with diameters of 8 nm have recently shown promising performance for biomedical applications. However, it is still unclear how the shape of cobalt clusters changes with size when reaching the nanoparticle range. In the present work, density functional theory calculations have been employed to compare the stabilities of two non-crystalline (icosahedron and decahedron) shapes, and three crystalline motifs (hcp, fcc, and bcc) for magic numbered cobalt clusters with up to 1500 atoms, based on the changes in the cohesive energies, coordination numbers, and nearest-neighbour distances arising from varying geometries. Obtained trends were extrapolated to a 104 size range, and an icosahedral shape was predicted for clusters up to 5500 atoms. Larger sized clusters adopt hcp stacking, in correspondence with the bulk phase. To explain the crystalline/non-crystalline crossovers, the contributions of the elastic strain density and twin boundary from the specimen surfaces to the cohesive energy of different motifs were evaluated. These results are expected to aid the design and synthesis of cobalt nanoparticles for applications ranging from catalysis to biomedical treatments.