Towards a morphology of cobalt nanoparticles: size and strain effects

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.


Introduction
In the last few decades, cobalt nanoparticles have been investigated and used mainly in catalysis [1,2] and magnetic data storage fields [3][4][5]. However, a long list of shortcomings of materials that have so far been implemented in magnetic hyperthermia treatments (e.g. iron oxides [6,7] and doped iron oxides [8,9]), including insufficient magnetic moments and too big diameters for the desired heating effects [10,11], combined with the high magnetisation of cobalt has opened up a completely new potential field of applications for this transition metal. Experimentally obtained heating rates for different materials [12] highlighted cobalt nanoparticles with diameters of ∼8 nm as one of the most promising candidates for novel cancer therapy. Unfortunately, compared to the mainly metallic oxide compounds considered previously for the same application, cobalt as a pure metal has a different drawback-it is highly reactive and easily oxidised even in mild conditions [13][14][15][16]. For that reason alone, additional modifications in the form of biocompatible surface coatings are necessary for its clinical implementation.
Applied coating could be organic or inorganic, in one or multiple layers, it could be linked to a drug carrier or even add advantageous plasmonic properties [17,18]. With so many possibilities, laboratory searches for optimum coatings are extremely time-consuming. Moreover, revealing the interactions between the base metal and coating material, and gaining an insight in the performance of such a complex system is a challenge in its own right. The increasingly predictive nature of computational investigations encourages the employment of computer modelling a priori, with experiments executed only Nanotechnology Nanotechnology 31 (2020) 195711 (15pp) https://doi.org/10.1088/1361-6528/ab6fe0 Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. for a few of the best combinations. However, for computer simulations of this type to be carried out, the exact geometries of the nanoparticles within desired size-range have to be known.
Understanding the structural, chemical, and physical properties of both small metal clusters and larger nanoparticles is obtained from a basic knowledge of their molecular and electronic structure. A first step is to determine the number of atoms in a cluster of a given size, their arrangement, and the final shape of the particle. This is followed by a systematic study of the dependence and evolution of the desired properties with the particle size and geometry to obtain the border transition sizes between different structural forms and to elucidate the gradual transformation from the finite molecular systems to bulk regimes.
Current experimental techniques allow the production of size-selected clusters with a narrow size distribution [19][20][21]. However, it is much more demanding to identify the atomic structure of generated clusters, which therefore remains an intriguing topic for both experimental and theoretical research. Thus far, computer simulations have found that three-dimensional isomers are always more stable than their planar counterparts for the smallest clusters of various metals [22][23][24]. Moreover, both theory and experiment have shown the preference of non-crystalline geometries (icosahedron and decahedron) over the crystalline shapes (fcc, hcp, or bcc) in clusters up to 50 or 100 atoms, with the crystalline/noncrystalline crossover point depending on the species. For palladium, platinum, and gold the transition from icosahedron to decahedron and fcc clusters was observed for N < 100 [25][26][27], but for rhodium and silver at N > 300 [28] and N > 400 [26], respectively. For an hcp ruthenium, the intersection of icosahedral and hcp structure was found to be at N=103 [29].
Despite the fact that for most transition and noble metals, non-crystalline to crystalline crossings were observed for relatively small cluster sizes, this is not the case for three ferromagnetic metals. The most favourable shapes for iron, nickel, and cobalt clusters with up to 800 atoms have been examined experimentally [30,31]. The geometry of iron clusters favours an icosahedral structure, which is, however, in competition with different crystalline motifs similar to the bulk phase and no clearly preferred shape has been confirmed. The exact geometries of smaller clusters of nickel and cobalt, up to 70 atoms, have also not been clearly identified; icosahedral features are less evident, and coexistence of isomeric species is presumptive. However, for larger clusters, a regular structure consistent with the filling of successive shells of atoms corresponding to an icosahedral cluster geometry was experimentally observed even for the clusters with 800 atoms. Moreover, in the case of nickel, an icosahedron was theoretically predicted to be the most stable structure up to 2300 atoms [32,33], which is amongst the largest size domains in which the icosahedral motif has been shown to be more stable than the bulk phase.
To date, theoretical studies on cobalt clusters are still incomplete and lack agreement. Early density functional theory (DFT) works were limited in terms of the cluster size and therefore only considered clusters with up to 10 or 15 atoms, mostly excluding geometry optimisation [34][35][36]. A later study performed on clusters up to 177 atoms, again without structural relaxation, took into account only fcc and bcc configurations [37], whereas tight-binding molecular dynamics studies have relied only on the experimental data of the cobalt dimer [38,39]. More recent works were carried out as spin-polarised relaxations without any symmetry constraints, but only for a few chosen sizes (2N 20 [40], N < 6 [41], and N=13, 55, 147 [42]). The most complete study so far has been conducted for medium clusters with up to 365 atoms in various structural motifs using first-principles electronic-structure calculations, but without considering decahedron or variations in fcc-shaped clusters [43]. Altogether, many cluster sizes have been disregarded, while for those discussed there is no consensus on the most stable geometries. Existing discrepancies among DFT calculations may be attributed to the choices of the basis sets and exchange-correlation functionals, and the incompleteness of the search for the global minimum by discarding some of the non-crystalline and/or crystalline motifs. Semi-empirical potentials and parametrised tight-binding methods drastically reduce the computational cost compared to DFT results, but with the loss of accuracy and reliability. The overall picture of the shape as a function of size is therefore still incomplete and the motif preferred by larger sized Co clusters remains an open question.
A suggested Co particle with 8 nm in diameter which yields the best heating power would constitute of more than 10 4 atoms and our main goal is, therefore, to identify the shape of clusters in the medium to large (N=10 2 -10 4 ) sizerange. Even the most advanced equipment struggles with clusters involving a few hundred atoms since they are too big to analyse by spectroscopic methods, and yet too small to analyse by diffraction methods [44]. On the other hand, a systematic theoretical study of larger cobalt clusters with up to several hundreds of atoms by accurate DFT calculations is still lacking, likely due to the hexagonal close-packed structure of the bulk material and costliness of calculations. Beyond a certain size, clusters are expected to exhibit bulkwise close-packed structures, despite the experimentally reported preference for icosahedral geometries up to 800 atoms. To the best of our knowledge, the crossover point between the icosahedral and hcp structural motifs has not as yet been determined.
In this work, DFT calculations have been carried out on cobalt clusters of different crystalline and non-crystalline structures containing up to 1500 atoms and results were extrapolated to capture behaviour of larger nanoparticles. The main interest of the study was to observe the size and strain effects caused by the appearance and intersection of different major surfaces and their impact on the stability of structural forms. Evolution of structural, bonding, and magnetic properties was also investigated. As we expect the transition from icosahedron to hcp motif to occur somewhere in the mediumsize range, with different strain contributions acting on the (in)stability of the non-crystalline structures, we have considered the following issues: (i) Can the chosen level of DFT provide enough accuracy for modelling cobalt clusters? (ii) Given the size (8 nm) and composition (monometallic cobalt) of the system of interest, what is the most stable cluster structure from an energetic point of view? and (iii) What are the effects of the cluster size on the overall contribution of the elastic strain energy and twin boundary from the specimen surfaces that characterise the cluster motifs?

Models and methods
Ideally, all cluster sizes should be tested, however, the number of possible structures and isomers grows non-linearly with the number of atoms, making such studies computationally prohibitively expensive. What can be done instead, and has already been reported in previous publications [45][46][47][48], is to consider a number of structures whose selection is motivated by energetic and morphological considerations. This choice is usually driven by the well-known recurrence of clusters which have a complete, regular outer geometry designated as full-shell clusters. These numbers are known as magic numbers, and they are given by the equations which differ for each shape [49]. In this study, multiple isomers (excluding bilayer structures) for the Co clusters ranging from 2 to 30 atoms (non-scalable regime) have been considered, while for the clusters with more than 30 but less than 1500 atoms (scalable regime), only the magic numbers of different motifs have been taken into the account. As the main focus is on the medium and large size range, a detailed search for low lying isomers of clusters under 30 atoms via global optimisation with a semi-empirical potential [50][51][52] was omitted, whereas up to ten different isomers were considered, depending on the number of atoms in a cluster.
Various levels of theory have been used when predicting the structures of transition metal clusters [23,[53][54][55][56][57][58]. Reliability of empirical potentials suffers from the neglection of the directional nature of d-d interactions and quantum effects arising from spin magnetism and orbital symmetry. For ab initio methods, structural predictions are challenging because of the significant contribution from localised d electrons and strong electron correlation. However, if the determination of both the geometry and electronic structure is to be achieved at the same level, which is desirable due to their strong interdependence, DFT proved to be a very efficient and reliable approach [42,[59][60][61][62][63][64]. Specifically for cobalt, it was shown that the lowest energy structures and properties of small clusters are not as sensitive to the level of the exchange-correlation functional used (including varying values of Hubbard and hybrid parameters), and results indicate that the differences between generalised gradient approximation (GGA), GGA+U, and hybrid methods are purely an electronic effect [65]. We have verified that the use of the GGA+U functional with proposed value of U=1.0 eV instead of GGA in the structural optimisation of small Co clusters would have had only a negligible effect in the form of a slight alternation in the relative energy difference of structural isomers, figure S1, supplementary information (SI) is available online at stacks.iop.org/NANO/ 31/195711/mmedia. While for small transition metal clusters with localised d orbitals GGA suffers from unphysical electron self-interaction, this deficiency diminishes as the cluster sizes approach bulk-like hybridisation and GGA is then more reliable for structural parameters and relative energetics. Therefore, while a comparative study of hybrid functionals versus GGA for absolute values of electronic properties for the small size region of Co clusters remains a future work, searches for the most stable structures of clusters with more than a few atoms using these high-level quantum chemistry methods would be very expensive. Because DFT-GGA is widely used, reported, and relied upon, it is worthwhile to find the structural trends at this affordable level, especially comparing the trends rather than concrete values. Thus, we consider that the choice of a functional and effective core potential used in current density functional calculations should be accurate enough to describe trends for cobalt clusters in the medium and large size range.
Spin-polarised DFT calculations were conducted on the clusters using the Vienna ab initio simulation package code [66] within the GGA exchange-correlation functional as developed by Perdew, Burke and Ernzerhof [67]. The core electrons up to and including the 3p levels of cobalt were kept frozen and their interaction with the valence electrons was described by the projector augmented wave method [68]. The kinetic energy cutoff of the wave function was set to 400 eV, and the Monkhorst-Pack k-point grid included only Γ-point. Relaxation of atoms was carried out without any constraints, with the convergence criteria for the change in the total energy and interatomic forces between two steps of 1.0×10 −6 (1.0×10 −3 for clusters exceeding 700 atoms) and 1.0×10 −2 eV Å −1 , respectively. The vacuum space in the unit cell was set to be more than 12 Å in each direction to avoid interactions with the periodic images.
In order to quantify the results, several quantities for the comparison of cluster stabilities have been introduced. The first is the cohesive energy, expressed as: with a to d being material and shape-dependent coefficients; strain enters along with the bulk-like cohesive energy in a, surface and twin boundary energies in b, edge and corner contributions in c and d, respectively; N stands for the number of atoms in the cluster. As a result of the competition between those contributions, different stable motifs occur in different size regimes. In terms of calculations, the cohesive energy per cohesive energy, roughly divided by the number of surface atoms, N 2/3 . Since the cohesive energy per atom in clusters/ nanoparticles is less negative than in the bulk, Δ(N) is positive and, in the limit  ¥ N , it tends to a constant value for crystalline nanoparticles, whereas it diverges for non-crystalline structures presenting volume strain. Stability of the small metal clusters is studied by plotting the second energy difference: 1, and N atoms, respectively. The peaks of the function represent relatively more stable clusters and correspond to the intensity anomalies or magic numbers of the mass spectrum.
The contributions of the elastic strain (for non-crystalline motifs) and twin boundary surface energies (for hcp and icosahedron) were evaluated based on the elastic constants of fcc cobalt (c 11 , c 12 , c 44 ), and twin boundary models [69], respectively. The twin boundary energy per unit area, g t , is expressed as: where A is the area of the twin boundary, and E original DFT and E twin DFT represent the calculated total energies of the original surface model and the twin boundary model, respectively.
The elastic strain densities for icosahedron and decahedron, W i and W d , were calculated by assuming a uniform elastic strain throughout the crystal, and elastic strain energies, E i and E d , expressed as follows: where r represents the length of the edge of icosahedron or decahedron [70][71][72][73][74][75]. More details on models used can be found in SI.

Small clusters-accessing the accuracy
The smallest clusters with 2 N 30 atoms provide an excellent test for the method for two reasons: first, the clusters have so few atoms that the number of (meta-)stable structures is still low enough for the global minimum to be easily accessible within a given approach. Second, since the planewave method has its foundations in considerations for infinite, extended systems, small clusters should be the most difficult to evaluate. Moreover, the biggest amount of data available covers clusters up to 30 atoms. Figure 1 shows optimised structures, average bond lengths, binding energies, average magnetic moments (per atom), and HOMO-LUMO energy gaps for the most stable shapes of clusters with 2-30 atoms. Expanded results for second isomers are provided in figure S2, SI. The obtained parameters are generally within a few % of the results of other parameter-free electronic-structure calculations (LDA [34], GGA [36,40], BLYP [41], DV-LSD [76]), with major differences only arising when considering the relative order of different isomers. Moreover, none of the parameterised methods, i.e. the tight-binding approach [39] or model potentials of Gupta [77], gives results that are systematically more accurate than the plane-wave approach. Furthermore, although there are discrepancies with other theoretical studies in terms of the most stable structures for some of the sizes, those predicted here are in good agreement with structures suggested in experimental works [78][79][80], with the 19-atom hcp structure being determined as more stable than the double icosahedron. For clusters ranging between 15 and 30 atoms, where sizes are very different from the values of closed-shell clusters with 13 or 55 (icosahedron, fcc)/57 (hcp) atoms, there is a strong competition between the two morphologies.
With the aim of identifying general properties of the clusters, their evolution with size, and accuracy of the results obtained, detailed descriptions of each isomer considered will not be provided (as most have already been treated elsewhere [39,40]). Instead of discussing the individual clusters, different quantities shall be introduced to reduce the collected information to a few key numbers. Figure 2(a) shows the comparison between calculated and experimental cohesive energies. Although DFT predicts a more uniform trend, the energy values are within the experimental range. To confirm this qualitative agreement, the second derivation of the energy as the stability function has been plotted in figure 2(b) against the number of atoms in the cluster. This function, which has maxima for particularly stable clusters, can be used to identify or, in this case, confirm the so-called magic numbers. The most pronounced peaks occur at N=6, 13, 19, which is in agreement with previously determined magic numbers for transition metals. An earlier study [40] did not capture the stability of the 13-atom cluster, which represents the first filled shell. Based on the occurrence of peaks for 6, 13, and 19 atom clusters in this study, it can be concluded that the cohesive energy values are qualitatively reasonable.
The magnetisation trend is shown in figure 2(c) as magnetic moment per Co atom for ferromagnetic ordering. Energies of clusters with antiferromagnetic ordering were calculated to be less favourable by, on average, 0.35 eV per Co atom. Experimental values provided are taken from the work of Knickelbein which investigated cobalt clusters in a similar size range (N=7-32) [82]. Underestimation of the magnetic moment as compared to the experiment could partially be caused by omitting the orbital moment, which can be rather large in clusters, in contrast to solids where it is usually strongly suppressed [83]. It was recently found that the orbital magnetic contributions that emerge in small cobalt clusters can even reach values that are an order of magnitude higher than those in the bulk [84]. Observable inconsistencies in the experimental and computed trend, however, can be improved for intermediate sizes (10 N 25) by introducing magnetic moments of the second most stable isomer, which confirms the possible coexistence of different isomers in the cluster beam that has been proposed experimentally [85]. Additionally, experimental measurements of magnetic moments are challenging at the smallest sizes: the cluster beam intensities decrease rapidly, thermalisation can be very difficult to achieve, and the net magnetisation of a cluster can change direction spontaneously due to thermal fluctuations. There are also temperature effects to be taken into account-a severe problem is still the determination of the cluster temperature and the distribution of the rotation frequencies, which are dependent on the experimental conditions [86,87]. As a result, variations in magnetisation values and trends can be found in the literature [88][89][90] making comparison rather difficult. Finally, magnetic moments of cobalt clusters approach the bulk value remarkably fast, even when compared to two of the metalʼs closest neighbours, iron and nickel [91].
Certain properties, such as the HOMO-LUMO gap values, are not experimentally available. Nonetheless, nickel clusters showed clear similarities to cobalt clusters in both shape and reactivity [78,92], making the comparison suitable. The trend of changes in the energy gap with the increase in the number of atoms in the cluster, shown in figure 2(d), resembles the trend for nickel clusters [93]. It is worth noting how, for both metals, gap values are significantly lower and clusters progress to metallic behaviour much faster than other transition metals [94,95], with the gap for cobalt clusters already closing at N=14−16.
Figures 2(e) and (f) summarise the quantities related to the structural properties, namely average coordination numbers and bond lengths. Two atoms are considered bonded if their interatomic distance is less than the average value of the first and second nearest-neighbour distances in the bulk (2.44 and 3.50 Å, respectively), which in this case is 3.00 Å. According to the coordination number, atoms are classified as inner if CN12, or surface atoms if CN < 12. Both properties oscillate mostly for the smallest clusters while increasing as a function of the size. Convergence towards the bulk values is obvious (2.44 Å nearest-neighbour distance and 12 coordination number), with the average coordination number saturating at a slower pace, due to most of the atoms still being characterised as surface atoms at such small sizes.
Taking everything into consideration, it was found that for the smallest most-critical clusters, the GGA results show satisfactory qualitative agreement with available theoretical and experimental information, making the method sufficiently accurate to evaluate larger sizes.

Large clusters-towards the behaviour of nanoparticles
Five kinds of motifs have been considered for clusters with up to 1500 atoms, two non-crystalline: icosahedron and decahedron, and three crystalline: hcp, fcc, and bcc, figure 3. By construction, the centre of the icosahedral and decahedral clusters is occupied by a single atom. Marks and Ino decahedra have been derived from the regular decahedron by truncating it by ten re-entrant (111) facets or by five (100) facets at the five twin boundaries, respectively. Chosen fcc structures are the atom-centred cuboctahedral and atom-and octahedron-centred truncated octahedral systems. Hcp clusters are truncated hexagonal bipyramids with cluster centres located at a single atom, a triangle (two possible combinations, with and without the central atom in the second layer), or interstitially between two (0001) layers. Two types of bcc clusters are represented: an atom-and an interstitial-centred clusters. Construction of differently centred motifs is represented in figure S3, SI, while structures of irregular decahedral geometries are shown in figure S4, SI.
The stabilities of the geometrical structures were compared based on their cohesive energies: the lower the cohesive energy, the more stable the shape. Figure 4(a) shows the cohesive energies for each motif as a function of N −1/3 , with the specific values of the magic numbers and accompanying energies given in table S1, SI. Dependencies are well represented by linear regression and extrapolated to large cluster sizes. For relatively small sizes (N < 100), the stability order is: icosahedron> truncated octahedron > hcp > cuboctahedron>decahedron> bcc. Differences in energies for the clusters with the same number of atoms between any two shapes are close to or less than 0.10 eV, except for crystalline bcc clusters which have stability lower by 0.40 eV compared to the most stable icosahedral morphology in the smallest size region. Decahedron and cuboctahedron are separated by minimal differences in energy; however, their stabilities do not intersect even with the extrapolation to larger diameters. Just under N=500, hcp motif surpasses fcc truncated octahedron in stability. After that crossover, the order remains unchanged with icosahedron being the most stable shape throughout the whole range of small and medium cluster sizes (N 5000), which is consistent with experimental findings of its appearance within clusters with up to 800 atoms. The icosahedron-to-hcp transition is predicted to happen at around N≈5500. An enlarged view of the intersection is represented in figure 4(b). Equations of the linear regression were used to get a crude approximation of the exact crossover point which is now suggested to be at N=5341. The intra-motif competition and stability of differently centred hcp, truncated octahedron, and bcc clusters, as well as of irregular decahedra is represented in figures S5 and S6, SI.
It should be noted that this type of fitting approach and a subsequent search for stability trends is constrained by assuming that the differences in stabilities of non-magic clusters of different structural motifs are close to the differences defined by the cohesive energies of magic clusters. By keeping in mind the experimental indications of the expected occurrence of more than a single geometrical shape at practically any size, sharp motif transitions as predicted by linearly interpolated energetics of magic numbered clusters should be taken only as a guideline to the predominating geometry. However, recent works that have sampled the energy landscape beyond the magic numbers found a qualitative agreement, where the most dominant cluster shape coincides with the one predicted by a simple 'magic number' fit [96][97][98]. The non-crystalline/crystalline distributions thus obtained for cobalt clusters with magic number of atoms represent a good reference point for defining crossover sizes between the structural motifs. Note that shape alternations could occur in reported stability windows, but they should nevertheless be expected to contain the highest proportion of the energetically most favourable structure. That being said, a spherical nanoparticle with diameter of 4.0 nm as cut from fcc Co bulk counts just under 8600 atoms (N=8586, figure S7, SI), and, with hcp morphology becoming successively more stable over icosahedron after the crossover point, by reaching 8 nm the divergence in the stability of the two would be significantly in favour of the hcp morphology and the proportion of hcp particles is therefore expected to far outweigh other motifs. However, confirmation through the continuous sampling of all cluster sizes remains a topic for future work.
To understand the shifts in the most stable structures with the increase in the cluster sizes, the effects of surface formation, twin boundaries, and elastic strain energy losses created at the expense of high surface-to-volume ratios have been estimated along with the evolution of structural (CN, Co-Co distances) and electronic (magnetisation) properties.
Each of the crystalline and non-crystalline motifs is characterised by specific Miller-index surfaces; of fcc shapes, cuboctahedron has six rectangular (100) and eight triangular (111) surfaces, truncated octahedron has six rectangular (100) and eight hexagonal (111) surfaces, while hcp bipyramids contain two (0001) and twelve (101̅ 1) hexagonal surfaces. Icosahedron and decahedron are built from twenty and ten (111) facets, respectively. The fcc (111) surface is wellknown as the surface with the lowest energy amongst low-Miller index surfaces of both hcp and fcc cobalt [99]. The stability order of four surfaces present in different cluster motifs is fcc (111) > hcp (0001) > hcp (101̅ 1) > fcc (100), (table S2, SI), which is connected to the coordination numbers, with atoms in the first two surfaces having CN=9, and CN=8 in the latter two. The dependence of average coordination number on the shape and size of cobalt clusters is shown in figure 5(a). Icosahedral clusters have the highest average CN for the whole size range as they contain only fcc (111) surface. Decahedral CNs are somewhat lower as a decahedron is built from half the number of (111) facets of icosahedron. Fcc (100) surface, characterised by CN=8, is introduced in cuboctahedron and truncated octahedron clusters. Since a cuboctahedron has a higher share of the lowercoordinated (100) surface atoms, the average CN is consequently lower when compared to both the non-crystalline shapes and the truncated octahedron. Hcp cobalt clusters have the lowest coordination numbers. Formation of an icosahedron is therefore favoured as fewer bonds need to be broken to achieve higher CN values, which lowers the energetical cost. The difference between the average CNs of all the motifs decreases with the growth of the clusters, eventually reaching the maximum of 12 bonds as in the bulk phase. Table 1 summarises the surface energies of magic-numbered cobalt clusters of each structure. The excess energy due to a free surface has two opposite effects: the energy cost related to the broken bonds and lowering of the coordination numbers, and the energy gain through the surface stress state. The surface energy of the hcp motif is the highest of the considered shapes for the similar sizes, which is directly connected to hcp cobalt having a higher cohesive energy (5.74 eV as calculated in this study) than the fcc bulk phase (5.44 eV as calculated in this study), causing increased energy losses when bonds are broken to form the hcp surfaces. This is also interconnected to variations in the energies of different hcp and fcc surfaces. As the cluster size decreases, the formation of surface facets has a significant effect on the overall stability, as the surface area dominates over volume. As a result, fcc motifs with mainly fcc (111) surfaces, i.e. truncated octahedron and icosahedron, become favoured. However, decahedral clusters retain significant surface areas in large clusters, with the 609-atom cluster still having more than 50% of atoms on the surface. Although this wide surface expanse has a comparatively lower surface energy, at the same time it demands that a high number of bonds are broken, which thus   Table 1. Surface areas, S (in %), and surface energies, γ (in eV Å −2 ) for magic numbered icosahedron, decahedron, cuboctahedron, truncated octahedron, and hcp cobalt clusters.

Icosahedron
Decahedron Cuboctahedron Truncated octahedron hcp makes it an unfavourable process as the size increases, and prevents better stability compared to cuboctahedron characterised by high-energy (100) facets. Average magnetic moments per atom for each shape are shown in figure 5(b) for ferromagnetic ordering. Because of the size of the clusters, different magnetic couplings (ferromagnetic, random and ordered antiferromagnetic, as well as a mixture of ferromagnetic and antiferromagnetic coupling within the cluster) were tested for cluster sizes between 105 atoms and 153 atoms. Ferromagnetic ordering was preferred by all motifs with rather large differences in energies of different couplings, e.g. for a 147-atom icosahedron, the ordered antiferromagnetic structure had an energy which was higher by 0.36 eV per Co atom when compared to the ferromagnetic one, while the mixed antifferomagnetic-ferromagnetic phase was less favourable by 0.30 eV per Co atom. For the hcp cluster with 153 Co atoms, the ordered antiferromagnetic structure was 0.46 eV less stable, and the mixed phase was 0.30 eV per Co atom less stable than the ferromagnetic phase. Initial random antifferomagnetic coupling was found to relax to the ferromagnetic ordering. Magnetisation behaviour with change in the size and shape reflects the relation between morphologies, coordination numbers, and surface shares. For sizes under 100 atoms, hcp crystalline motif has the lowest average magnetic moment despite its low average CN owing to the natural bulk packing. In medium size region, decahedron owes its high magnetic moment to the wide surface areas, while magnetisation of cuboctahedron arises due to the open geometry of fcc (100) surface, which, alongside low CNs, showcases moments of 1.82 m B for the first layer atoms (tables S2, SI). A detailed insight into the electronic structure differences that govern the magnetic behaviour of icosahedral and truncated octahedral Co particles can be found in a recent study [100]. Hcp shaped clusters approach the bulk magnetisation value with the increase in the size.
In addition to obvious differences in the energetics of appearing surfaces, there are additional factors induced by specific surface geometries and intersections that may be responsible for shape crossovers. The fcc (111) surface may be characterised by the lowest energy and the highest CN, but tetrahedral arrangements in the icosahedron and decahedron, while minimising the surface-to-volume ratio, at the same time introduce the elastic strain contribution. The elastic strain energy density, W, is expressed in terms of elastic constants c 11 , c 12 , and c 44 of fcc cobalt (tables S3, SI), and is calculated to be 1.95×10 −4 eV Å −3 for decahedron and 2.92×10 −3 eV Å −3 for icosahedron. Despite the known deficiencies of DFT when computing elastic properties of materials and variations in existing experimental data, the obtained elastic strain energy densities have the same order of magnitude and values that are close to the those published previously (3.00 × 10 −4 eV Å −3 for decahedron and 1.75 × 10 −3 eV Å −3 for icosahedron [69]). The overall contributions to the cohesive energies of the decahedral and icosahedral cobalt clusters are expressed in terms of elastic strain energy in table 2, and have the order of 10 −3 and 10 −2 eV per atom, respectively, with the values increasing with the growth in the size of the cluster. For both icosahedron and decahedron, and for all the sizes, elastic strain energies are positive, meaning they destabilise the noncrystalline clusters.
Another consequence of small volume packing that can either stabilise or destabilise a cluster, depending on its shape, is the presence of twin boundaries at surface intersections. To estimate the twin boundary energy per unit area for different motifs, twin boundary models of fcc (111) and hcp (101̅ 1) surfaces have been used, figure 6, with the calculated twin boundary surface energies collected in table 3. The obtained energies per unit area for hcp (101̅ 1) twin boundaries have positive values, indicating that this behaviour is undesirable and slightly destabilises hcp shapes-similar to what has been observed before for gold [101]. In contrast, for the fcc (111) surface the twin boundary energy per unit area is negative, therefore stabilising icosahedral and decahedral motifs. This behaviour is the opposite to what has been seen for transition metals which are naturally present in the fcc bulk phase [102], due to the similarity of the local structure near the fcc (111) twin boundary to the hcp (0001) stacking, which is the native bulk arrangement for cobalt. The twin boundary energy remains relatively unchanged with increasing periodic units, suggesting that its influence will eventually diminish as the cluster size increases.
Finally, volumetric contributions as a consequence of dense cluster packing and exposure of wide surface areas to vacuum also play a significant role in shape stability. Figure 7 shows the average interatomic distances between the nearestneighbour, surface-surface, surface-inner, and inner-inner Co atoms as a function of -N 1 3 for the motifs considered. Generally, the overall average interatomic distance is less than the bulk 2.44 Å and it increases as the cluster size grows, relatively quickly approaching the bulk value for all motifs. Distances between the inner cobalt atoms are the only ones enhanced compared to ideal 2.44 Å for all fcc and hcp crystalline shapes, but also for decahedron, while distances that include one and two surface atoms are up to 3.8% and 1.5% shorter than the bulk value, respectively.
Icosahedra of all sizes, on the other hand, have surfacesurface distances which are elongated more than the innerinner distances, with the elongation of the latter appearing only after the cluster size has reached 500 atoms. On contrary, the only 'inconsistency' from crystalline geometries in the trend of decahedron is for the smallest 105-atom cluster, where the distances between surface and inner Co atoms exceed that of the two surface Co atoms. Overall, for all shapes the variation in the distances between the inner cobalt atoms is minor compared to the distances which include at least one surface cobalt atom. Thus, the interatomic distances of cobalt atoms in the surface layer affect the average interatomic distance, and hence the stability of the motifs, the most. This is the most obvious for icosahedral clusters which are limited by (111) close-packed facets only, and the surface energy gain is optimised at the expense of a volume contribution. Consequently, inter-shell bonds are substantially compressed, and intra-shell bonds are expanded, imposing high strain for icosahedron structure, which confirms indications that it can only be favoured at sizes where the number of atoms cannot generate significant volume stress. Absence of any further crossovers between motifs at large cluster sizes is now explained through both crystalline and decahedral geometries obeying distance constraints to minimise the volumetric losses.
Taking everything into consideration, the highest coordination number and the lowest surface energy together with the negative twin boundary energy make icosahedral clusters the most stable for a wide range of cluster sizes. However, these contributions become less effective with increasing size, when Co-Co atom distances in icosahedron start to substantially deviate from the average values of the bulk. Moreover, the surface area naturally decreases as the clusters grow, reducing the contributions of surface energies in all motifs, as shown in figure 8. While the extent of the decrease in the surface energy of both crystalline and non-crystalline fcc shapes also reduces with cluster growth, the hcp surface energy drops quicker as the clusters approach medium and large size-range. Additionally, icosahedron and decahedron seem to be reaching a plateau already at N ≈ 1000, contrary to crystalline shapes. Nevertheless, even for clusters with more than 1000 atoms, there is still a noticeable difference in the surface energies of hcp and icosahedron motifs, which leads to a conclusion that volume and elastic strain become dominant factors in the (in)stability of the icosahedron only at very large cluster sizes, where the surface energies of the two shapes become less distinct.
Compared to another hcp metal, ruthenium, where the stability of hcp clusters already overcomes the icosahedral shape at N=103, cobalt requires significantly larger sizes before it experiences the same change. Figures 9(a) and (b) show the relation between size and strain effects of the two hcp metals. In the case of ruthenium, the surface energies of icosahedron and hcp clusters are close to each other even for sizes under 100 atoms. At the same time, the surface energies of cobalt icosahedral clusters are lower and of hcp higher, respectively, than those of the icosahedral and hcp ruthenium clusters. Combined with the higher values of the elastic strain for Ru icosahedron clusters, the surface contribution is, in contrast to cobalt, expected to become negligible at much lower sizes, as shown with a significant difference in the crossover point.

Conclusions
The dependence of the stability and properties of nanoparticles on their shape paves the way to the understanding of the intrinsic link between their versatile functionalities and morphology. Various motifs of cobalt nanoparticles, including crystalline and non-crystalline, atom-, triangle-, and interstitial-centred configurations have been modelled to reveal the transition sizes between different structural forms and the surface and energy strain contributions to the (in) stability of each of the shapes. DFT calculations within GGA approach have shown to be sufficiently accurate in predicting the stability of cobalt clusters and they were employed to define the size and strain  effects that govern the formation of different shapes as the number of atoms increases. The icosahedron was determined as the predominant motif for clusters up to 5500 atoms by virtue of the low surface and negative twin boundary energies. Its high average coordination number also contributes significantly, with volumetric and elastic strain effects becoming dominant only at larger sizes. Once the surface energies of icosahedron and hcp clusters reach comparable values, deviation of the interatomic distances from expected bulk values imposes a high strain to icosahedral clusters, thereby making the natural hcp arrangement preferred.
In conclusion, for larger clusters the minimisation of the internal strain energy is the most important factor. The strain energy is at its lowest in hcp clusters due to the bulk crystal nature, while it becomes the main destabilising factor for noncrystalline shapes. Considering the decahedral motif, the energetic cost of having larger surface areas is too great unless exceptional electronic configurations exist, which, for cobalt, is not the case.