Structure, thermodynamics, and rearrangement mechanisms in gold clusters—insights from the energy landscapes framework

Abstract

We consider finite-size and temperature effects on the structure of model Au_N clusters (30 ≤ N ≤ 147) bound by the Gupta potential. Equilibrium behaviour is examined in the harmonic superposition approximation, and the size-dependent melting temperature is also bracketed using molecular dynamics simulations. We identify structural transitions between distinctly different morphologies, characterised by various defect features. Reentrant behaviour and trends with respect to cluster size and temperature are discussed in detail. For N = 55, 85, and 147 we visualise the topography of the underlying potential energy landscape using disconnectivity graphs, colour-coded by the cluster morphology; and we use discrete path sampling to characterise the rearrangement mechanisms between competing structures separated by high energy barriers (up to 1 eV). The fastest transition pathways generally involve metastable states with multiple fivefold disclinations and/or a high degree of amorphisation, indicative of melting. For N = 55 we find that reoptimising low-lying minima using density functional theory (DFT) alters their energetic ordering and produces a new putative global minimum at the DFT level; however, the equilibrium structure predicted by the Gupta potential at room temperature is consistent with previous experiments.

---

I. Introduction

Gold has been one of the most important elements in cluster science, providing a model system for exploring fundamental questions,^1,2 and offering a range of size- and structure-dependent properties useful in various applications.³ Some novel properties, such as the unexpected catalytic activity,⁴ have been linked to specific cluster morphologies, and the perceived structure–property relationship has motivated many studies^5–8 of the atomic-level structure evolution. However, the diversity of morphologies and quasimelting^9,10 observed in gold nanoparticles has hindered systematic understanding, especially in discerning the equilibrium picture and characterising the rearrangement mechanisms at finite temperatures. In the present contribution we use the energy landscapes framework¹¹ to shed new light on the polymorphism in model gold clusters.

One of the first and still widely studied gold nanoparticles is the “Schmid Au₅₅” cluster,^12,13 formulated as Au₅₅ (PPh₃)₁₂ Cl₆. The Au₅₅ core geometry was first characterised as cuboctahedral,^13,14 illustrated in Fig. 1a, with single-crystal face-centred cubic (fcc) atomic ordering. However, this initial interpretation was criticised by Vogel et al.,¹⁵ who found that a model with icosahedral Au₅₅ core produced a better fit to the available X-ray powder diffraction data. More recently, Pei et al.¹⁶ used density functional theory to show that many quasi-icosahedral, decahedral, and disordered core structures are energetically favoured over the closed-shell cuboctahedron. To the best of our knowledge, the debate over the Au₅₅ core geometry has still not reached a definitive resolution, exacerbating interpretation of the unusual oxidation resistance¹⁷ and high catalytic activity¹⁸ of the naked Au₅₅ cluster derived from Au₅₅ (PPh₃)₁₂ Cl₆. The intriguing chemistry of the naked Au₅₅ may well be at least in part due to a particular geometry,¹⁷ but the geometry is unlikely to be cuboctahedral, because previous theoretical calculations^19,20 show that the expected lifetime of bare Au₅₅ cuboctahedra is too short to be observed under an electron microscope, even in the presence of a substrate.¹⁹

Fig. 1 The structure of some Au₅₅ isomers from Table 1: (a) Cub; (b) 5; (c) 8; (d) 1.

In another theoretical contribution, Garzón et al.^21,22 found several amorphous low-energy structures for the naked Au₅₅, with “amorphous” signifying distorted icosahedral order (see Fig. 1b and c), which can exhibit a degree of chirality²³ and enantioselective properties.²⁴ The lowest-lying amorphous isomer (Fig. 1b) was proposed as the global minimum (GM). Indeed, distorted²⁵ or amorphized¹⁰ icosahedra and the partially chiral⁸ Au₅₅ structures have been observed under an electron microscope. However, using the same empirical model for Au₅₅ as Garzón et al.,²¹ Bao et al.²⁶ found a lower energy structure with fcc order, shown in Fig. 1d, which had previously²⁷ been identified for a different interatomic potential. These theoretical predictions suggest that microscopy studies may not be accessing the lowest energy structure, perhaps due to finite temperature effects, an underlying substrate,²⁸ the electron beam, or some other factors. Recently,²⁹ we performed global optimisation of Au₅₅ for a different parametrisation³⁰ of the model used by Garzón et al.²¹ and Bao et al.,²⁶ and again found the GM to be fcc, with amorphous isomers dominating at finite temperatures, suggesting that thermal effects are indeed contributing to the discrepancy. In the present study we take a closer look and find that the mean occupation probability of the chiral isomer observed by Wang and Palmer⁸ is actually comparable to the equilibrium occupation probability of the model bound by the same potential as in ref. 21 and 26.

Prediction of an fcc GM for Au₅₅ by different empirical models has not yet been adequately tested at higher levels of theory. While DFT has been used to show that the lowest-lying amorphous²² and fcc³¹ isomers are individually more stable than other structures, as far as we know the two isomers have not been compared directly using the same DFT functional. Hence, the question of what is the true ground-state morphology of naked Au₅₅ still remains unresolved from a theoretical viewpoint. We address this issue through a more comprehensive exploration of structures and direct comparison between different levels of theory, producing a new putative GM for Au₅₅ at DFT level. We also verify that high-symmetry morphologies (such as the Mackay icosahedron,³² Ino decahedron,³³ and cuboctahedron) with closed geometric shells are not as stable as one might expect. Depending on the level of theory used, this destabilisation of high-symmetry structures has been linked to either the range of the interatomic potential^21,34,35 or relativistic effects.^31,35

Note that identifying the GM is necessary for an accurate description of the equilibrium behaviour within a given model, which is why global optimisation of gold (and other metal) clusters remains an active area of research.³⁶ However, the GM alone is not sufficient to explain the finite-temperature behaviour and morphological changes observed in experiments^5,7,8 and simulations.^29,37 In fact, finite-system analogues of solid–solid phase transitions have been reported for many metals,^38,39 sometimes well below the size-dependent¹ melting temperature, but systematic theoretical analysis of this phenomenon for a range of cluster sizes has been performed only for Lennard-Jonesium.^40,41 This omission is partly due to technical difficulties, because the relatively long time scales associated with morphological rearrangements below the melting temperature range cannot be easily accessed using conventional simulation methods. The energy landscapes¹¹ framework, on the other hand, provides a powerful approach to studying such rare events, complementing more conventional methods. This framework combines a variety of optimisation and search techniques, statistical mechanics, unimolecular rate theory, and often exploits the harmonic approximation to describe the global thermodynamics and kinetics of complex systems such as atomic clusters.

In the present contribution we explore the potential energy landscape of model Au_N clusters (30 ≤ N ≤ 147), focusing on their equilibrium thermodyamics in the harmonic superposition approximation,⁴² which accounts for configurational and vibrational entropy. We identify a number of solid–solid transitions in morphology, which arise from the competition between close-packed (single-crystal fcc or lamellar twinned), decahedral, and distorted icosahedral motifs. In selected representative cases we also find the fastest transition pathway between competing motifs. These pathways rarely exhibit the highly cooperative rearrangements identified in other systems,^43,44 but rather involve more localised distortions and are reminiscent of the melt-freeze scenario described by Koga et al.⁷

The outline of this paper is as follows: in section II we define the Gupta potential, describe the relevant methods from the energy landscape framework (including molecular dynamics simulation), detail our approach to classifying atomic-level structure, and give the methodological details of our DFT calculations. All the results are discussed in detail in section III, where we first focus on the structure and thermodynamics of the naked Au₅₅ cluster, then examine other cluster sizes, and in the end elucidate some rearrangement mechanisms for selected cases. A summary of key findings and conclusions is given in section IV. For completeness, in ESI† we provide: (i) a database of local minima on the Gupta potential energy landscape for all the Au_N clusters considered, including the coordinates of the putative global minima; (ii) kinetic transition networks for N = 55, 85, and 147; and (iii) the coordinates for the low-lying Au₅₅ structures reoptimised using DFT.

II. Model and methods

We model the binding in Au_N clusters of size N using the Gupta potential^30,45,46where r_ij is the distance between atoms i and j, and the remaining parameters are ξ = 1.855 eV, A = 0.2197 eV, p = 10.53, q = 4.30, and r₀ = 2.88 Å, as in ref. 35. Note that (1) is based on the second moment tight binding approximation, which remains a standard choice for atomistic modelling of metal nanoparticles,⁴⁷ and in the present study the potential is not truncated (hence no cutoff radius).

A. Harmonic superposition approximation

To explore the underlying potential energy landscape for a given N, we first perform basin-hopping^48,49 global optimisation aided by systematic surface refinement,⁵⁰ with atom-vacancy swap candidates scanned in parallel (as implemented in GMIN⁵¹). More than 10⁴ lowest-lying minima accumulated during global optimisation provide the input for a harmonic superposition analysis of the thermodynamics.⁴² This analysis involves writing the canonical partition function aswhere the sum runs over the local minima, E_m is the potential energy of minimum m, and k_B is the Boltzmann constant. The degeneracy factor^11,29subsumes entropic contributions, represented by the point group order (o_m) and the geometric mean normal mode vibrational frequency ([small nu, Greek, macron]_m). Here, κ = 3N − 6 is the number of vibrational degrees of freedom, h is the Planck constant, andwhere λ_i^(m) are the positive, non-zero eigenvalues of the mass-weighted Hessian (the dynamical matrix) for minimum m. The occupation probability of each minimum is given by

p_m(T) = g_me^−E_m/k_BT/Z(T), (5)

and the vibrational contribution to the heat capacity iswhere , and the constant kinetic contribution has been omitted. Since the p_m are additive, we can define collective occupation probabilities , where is a subset of minima.

B. Kinetic transition networks

For selected cluster sizes we construct a kinetic transition network, whose nodes represent local minima on the potential energy landscape, with edges connecting pairs of adjacent minima separated by a single intervening transition state. Pairs of low-lying minima obtained from global optimisation are first systematically connected using a Dijkstra-based approach⁵² implemented in OPTIM.⁵³ The pairwise connection attempts are performed with the doubly-nudged⁵⁴ elastic band algorithm,⁵⁵ with the intervening transition state candidates tightly converged using hybrid eigenvector following,⁵⁶ and intervening local minima converged using the limited-memory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) algorithm of Liu and Nocedal.⁵⁷ The resulting database of stationary points is further extended using discrete path sampling,⁵⁸ exploiting previously described⁵⁹ strategies implemented in PATHSAMPLE.⁶⁰ We visualise these networks in the form of disconnectivity graphs,⁶¹ which provide a revealing picture of the landscape topography.

To shed light on the rearrangement mechanisms between competing structures, we consider the possible connecting pathways in the corresponding network and identify the one with the largest contribution to the steady-state rate constant.^58,59 This pathway, referred to as the “fastest” for a given temperature, corresponds to the lowest-energy pathway in the limit of zero temperature.

The kinetic lifetime of individual minima is calculated using harmonic transition state theory.⁶² That is, the escape rate from a minimum m via an adjacent transition state s is

where the pre-exponential frequency factor is given bywith i and j spanning the positive (non-zero) eigenvalues of the mass-weighted Hessian matrix evaluated at m and s, respectively. The temperature-dependent lifetime of m is then defined aswhere the sum is over all the directly connected transition states. Note that ν₀^m→s can be smaller than [small nu, Greek, macron]_m defined in (4), which means that above a certain temperature the lifetime τ_m will be shorter than the geometric mean period of vibration, signalling a likely breakdown of transition state theory.

C. Atomic-level structure classification

The atomic-level structure of local minima is characterised using common-neighbour analysis (CNA),⁶³ with the nearest neighbours defined by a cut-off distance of 3.5 Å. First, the local order around each atom is classified (using the scheme of Hendy and Doye⁶⁴) as either icosahedral (ico), hexagonal closed-packed (hcp), or face-centred cubic (fcc), with the latter class including well-defined (100) and (111) facets, but not other surface features such as islands, re-entrant grooves, etc. Unclassified atoms are simply labelled as ambiguous (amb). We then characterise the overall cluster morphology, i.e. the motif, as either icosahedral (ICO), decahedral (DEC), face-centred cubic (FCC), twinned (TWI), hexagonal close-packed (HCP), or otherwise ambiguous (AMB), using a sequence of simple criteria similar to those in ref. 28. First, if an ico atom is neighbour to more than six other ico atoms, then it is regarded as an icosahedral centre and the motif is classified as ICO. Otherwise, if every ico atom is bonded to at most two other ico atoms, and if the total ico atom count is one less than the number of bonds with CNA signature 555, then there is one (local) fivefold symmetry axis (i.e. a decahedral spine) and the motif is labelled as DEC. Otherwise, if the number of ico atoms is zero and the combined number of hcp and fcc atoms exceeds the number of amb atoms, then the motif is either FCC if the number of hcp atoms is zero, HCP if the number of fcc atoms is zero, or else it is labelled TWI for “twinned”, though we make no attempt to distinguish between twin planes and stacking faults. Structures failing to satisfy all these criteria are classified as AMB, which is quite a broad category, often capturing minima that are similar to some from the other (more precisely defined) motifs.

D. Molecular dynamics simulations

To complement our harmonic superposition analysis, the thermal behaviour of selected isomers is simulated using canonical molecular dynamics with a velocity Verlet scheme coupled to a Langevin thermostat.⁶⁵ We used a damping coefficient of 5 ps⁻¹ (equivalent to a time constant of 20 fs), and the Verlet integration time step was set to 10 fs. Each simulation ran for 2 × 10⁷ time steps or more, with the first 4 × 10⁶ steps treated as equilibration and discarded. Over the remaining 160 ns we computed the time- and atom-averaged Lindemann index⁶⁶where 〈…〉_t indicates a time average of the quantity within the angle brackets. Decreasing the integration time step had no significant effect on δ.

E. Density functional theory calculations

Selected low-lying Au₅₅ minima were reoptimised using the Quantum Espresso package⁶⁷—a density functional theory based plane wave code. The exchange–correlation potential was described self-consistently within the generalized gradient approximation (GGA) throughout the Perdew–Burke–Ernzerhof (PBE) functional.⁶⁸ The Rabe–Rappe–Kaxiras–Joannopoulos ultrasoft pseudopotential was used to model the valence electron-nuclei interactions. The energy cut-off for the plane-wave basis set was 40 Ry with a charge density cut-off of 360 Ry. The Au electronic configuration considered was 5d¹⁰6s¹. All the calculations were performed at the Gamma point only in a cubic simulation box of at least 29 Å (deemed sufficiently large). Electronic ground state optimisation was performed starting from the ionic structure relaxed at the empirical (Gupta) level. A Gaussian smearing was used with effective electronic temperature of 27.2 meV. The spin–orbit coupling was not included, because it is not expected to have a significant effect on the energy and geometry of Au₅₅ minima.⁶⁹

III. Results and discussion

A. The naked “Schmid Au₅₅”

The disconnectivity graph for 500 lowest-lying minima of the Au₅₅ cluster is shown in Fig. 2. Recall that the vertical axis corresponds to the potential energy, each branch ends at a local minimum, and each node joins minima that can be interconverted without exceeding the energy of the node (thus providing an upper bound on the lowest-energy barrier between them). The energy of the nodes (but not branch endpoints) has been discretised in regular steps for clarity, revealing multiple competing funnels: one narrow funnel comprised of minima with close-packed structure (i.e. FCC and TWI motifs), and the others dominated by the AMB motif. Isolated minima of DEC and ICO motifs are also present, and the colour-coded structure of the lowest-lying minimum for each motif is illustrated below the disconnectivity graph. Note that many low-lying ICO and AMB minima resemble icosahedra with discernible triple rosette-like⁷⁰ defects and other surface distortions. However, single and double rosettes are still not particularly favourable, though we expect them to be stabilised by a smaller impurity atom.²⁹

Fig. 2 Colour-coded disconnectivity graph for the 500 lowest-lying minima of the Au₅₅ cluster. Branches leading to minima of icosahedral (ICO) motif are in red, decahedral (DEC) in green, face-centred cubic (FCC) in blue, twinned face-centred cubic (TWI) in magenta, and the remaining ambiguous (AMB) morphologies in black. The node corresponding to isomer 6 in Table 1 is marked by a square, and isomers 7 and 8 are circled. Ball-and-stick representation of the lowest-lying minimum for each motif is shown from two angles, with atoms colour-coded by the local environment: face-centred cubic (fcc) in gold, hexagonal close-packed (hcp) in white, icosahedral (ico) in red, and ambiguous (amb) in grey. The stick “bonds” are defined by a cut-off distance of 3.5 Å.

More details for the sixteen lowest-lying minima are given in Table 1, confirming the GM found by Bao et al.,²⁶ and showing that the lowest-lying AMB minimum found by Garzón et al.^21,22 is (at best) fifth lowest overall in the given model. The eighth isomer is identified as the chiral structure imaged by Wang and Palmer.⁸ The lowest-lying DEC minimum is fourteenth overall, with the fivefold disclination significantly off centre—in contrast to the Ino³³ and Marks⁷¹ decahedra, but reminiscent of some pentagonally twinned structures reported for lead nanoparticles.⁷²

Table 1 Potential energy relative to the GM (i.e. ΔE = E − E^GM), motif, point group (PG), the geometric mean normal mode vibrational frequencies (v̄) in terahertz, and the number of distinct adjacent transition states (N_t.s.) for the sixteen lowest-lying minima and three high-symmetry structures of Gupta Au₅₅. Each minimum's occupation probability (p_m) and lifetime (τ_m) in seconds is calculated at room temperature (k_BT = 26 meV) using a database of more than 3 × 10⁵ minima and 4 × 10⁵ transition states. The relative energy (ΔE_DFT) of each isomer reoptimised at the DFT level is also given, and the chiral isomer (ranked eighth) imaged by Wang and Palmer⁸ is marked by a dagger

	ΔEG (eV)	Motif	PG	v̄ (THz)	N t.s.	p m	τ m (s)	ΔEDFT (eV)
1	0	FCC	C 1	2.16602	4295	0.0026	2 × 10^−12	0
2	0.021294	FCC	C s	2.16593	435	0.0006	3 × 10^−11	−0.234
3	0.027287	TWI	C 1	2.16607	202	0.0009	8 × 10^−12	−0.361
4	0.030151	TWI	C 1	2.16641	160	0.0008	1 × 10^−11	−0.439
5	0.035735	AMB	C 1	2.09791	664	0.1050	2 × 10^−12	−0.503
6	0.037537	TWI	C 3v	2.16414	230	0.0001	9 × 10^−11	−0.039
7	0.039020	AMB	C 1	2.09591	276	0.1086	1 × 10^−11	−0.391
8^†	0.039037	AMB	C 1	2.09336	555	0.1313	1 × 10^−12	−0.551
9	0.043134	AMB	C s	2.13989	259	0.0017	5 × 10^−11	−0.759
10	0.055655	AMB	C 1	2.09451	176	0.0632	1 × 10^−12	−0.544
11	0.056596	TWI	C s	2.16562	198	0.0002	3 × 10^−11	−0.729
12	0.066701	TWI	C s	2.16574	86	0.0001	2 × 10^−12	−0.759
13	0.069885	AMB	C 1	2.12300	125	0.0044	8 × 10^−14	0.210
14	0.073645	DEC	C s	2.15466	339	0.0002	2 × 10^−11	0.449
15	0.088317	AMB	C 1	2.10861	164	0.0062	9 × 10^−14	−0.500
16	0.093786	ICO	C 1	2.07067	625	0.0899	2 × 10^−11	0.030
Mac	0.642895	ICO	I h	1.98361	46	0.0000	1 × 10^−10	0.644
Ino	1.070310	DEC	D 5h	2.10213	31	0.0000	4 × 10^−14	1.084
Cub	1.194899	FCC	O h	2.10972	29	0.0000	5 × 10^−19	1.748

---

While the GM structure of Au₅₅ is fcc with point group C₁, the symmetric cuboctahedron is significantly higher in energy (by 1.2 eV). This disparity can be explained by differences in surface packing: the 42 surface atoms in the cuboctahedron form only (100) facets, which are known³⁵ to be particularly unfavourable in the present model for gold; but in the GM structure the 45 surface atoms are more close-packed and exhibit mainly (111) character. Also, in agreement with previous^21,37 reports of distorted icosahedral order in Au₅₅, snapshots of the AMB minima in Fig. 2 exhibit triangular close-packed facets and multiple fivefold disclinations, which outline the tetrahedral units expected in Mackay³² icosahedra. The ideal Mackay icosahedron, on the other hand, is not even in the lowest-lying 10⁵ minima (in a database of more around 3 × 10⁵ minima), though it is the most favourable among the closed-shell high-symmetry shapes.

The disconnectivity graphs in Fig. 2 and 3 (discussed below) help us to visualise the landscape topography and to identify the funnels associated with competing motifs, but these graphs do not really show how the competition between different morphologies is manifested under thermal conditions. To examine finite-temperature effects we consider the heat capacity C_V and the collective occupation probability P^XXX_occ of each motif (XXX) as a function of k_BT > 0, plotted in Fig. 4.

Fig. 3 Disconnectivity graphs for 10³ minima of Au₈₅ (top) and 10⁴ minima of Au₁₄₇ (bottom) with ball-and-stick representations of the lowest-lying minimum for each competing motif. The colour-coding and nomenclature are same as in Fig. 2, and only a subset of atoms is highlighted for clarity.

Fig. 4 Motif occupation probabilities P_occ (stacked on top of each other), vibrational heat capacity C_V, and time-averaged Lindemann index δ plotted against temperature for Au₅₅ (top), Au₈₅ (middle), and Au₁₄₇ (bottom). Five motif shares (FCC, AMB, ICO, TWI and DEC) of P_occ are colour-coded as in Fig. 2 and 3. The chiral (eighth in Table 1) Au₅₅ isomer's occupancy is highlighted in grey, splitting the black segment in two, since contributions from individual minima of a given motif are stacked in the order of increasing potential energy. Net C_V is represented by a thick black line, while thinner black lines correspond to subsets of 10, 10², 10³, 10⁴, and 10⁵ lowest-lying minima in our database (sorted by potential energy). The Lindemann index δ was calculated from molecular dynamics simulations at temperatures marked by the datapoints, and datasets with the same starting configuration are colour-coded by the initial motif and traced by a dashed line to guide the eye.

For Au₅₅, P^FCC_occ is the highest below 70 K (k_BT < 6 meV). This high occupancy is associated mainly with the GM structure, while other low-lying FCC and TWI isomers combined have near-zero occupation probability at all temperatures, largely due to their low vibrational entropy. P^AMB_occ rapidly grows for k_BT in the range 6 ± 1 meV (T ≈ 70 K), with several AMB minima reaching comparably high occupation probability, which is illustrated by singling out the contribution from the experimentally observed⁸ chiral isomer (eighth in Table 1 and visualised in Fig. 1c). Recall that a similar transition at a slightly higher temperature (about 90 K) has been identified for a different set of model parameters.²⁹ In both cases, when P^AMB_occ supplants P^FCC_occ as the highest value, the crossover temperature coincides with a well-defined peak in the heat capacity. Hence, we characterise the finite-system analogue of a solid–solid like phase transition, where the low-energy phase is represented by a single local minimum of FCC motif, and the high-energy phase comprises multiple AMB isomers. Also note that the absence of magenta and green for Au₅₅ in Fig. 4 is due to very low occupancy of DEC and TWI motifs for the entire temperature range considered.

Interestingly, the room-temperature occupation probability of the chiral AMB isomer is 0.13, which is comparable to the occurrence frequency of about 0.1 inferred from experimental data of Wang and Palmer.⁸ This agreement between theory and experiment can be rationalised by taking into account the lifetime of individual isomers (see Table 1). The estimated lifetimes are below a typical vibrational period at room temperature, indicating interconversion among multiple isomers on a timescale significantly shorter than the experimental imaging time. Indeed, Wang and Palmer⁸ acknowledge that their images could be a superimposition of multiple isomers, which would explain why the average occurrence frequency of the chiral isomer observed under the microscope is comparable to the equilibrium occupation probability in Table 1. This rapid fluctuation between several different isomers could also be interpreted as quasi-melting.⁹

Fig. 4 shows that the ICO motif hardly features in Au₅₅ at low temperatures, but its occupation probability steadily increases over a temperature range that roughly coincides with a second and more dominant peak in the heat capacity. The time- and atom-averaged Lindemann index δ also increases from below 0.1 to above 0.3 in that temperature range, indicating a finite-system analogue of a solid–liquid phase transition (i.e. melting). The agreement between MD results and the harmonic superposition approximation is encouraging, with the latter formulation revealing interesting changes in the atomic-level structure of the molten Au₅₅. The rise of P^ICO_occ with k_BT in the range 30 meV < k_BT < 80 meV, with P^AMB_occ decreasing yet remaining significant, suggests a gradually growing preference for local (poly)icosahedral order (i.e. increasing fivefold-disclination density) in the melted region.

To cross-check our analysis of the Au₅₅ cluster, we reoptimised the structure of Gupta minima in Table 1 using DFT. The relaxation produces fairly minor geometric changes, mainly via uniform expansion or compression of the isomers. However, there is significant re-ordering of the energies: the fcc GM predicted by Gupta is not even in the ten lowest at the DFT level, where the new putative GM of point group C_s is obtained by relaxing isomers nine and twelve (in Table 1). The DFT GM structure is shown in Fig. 5, illustrating its amorphous nature and revealing two voids in the subsurface region (see Fig. 5b), consistent with the known propensity of gold clusters to form cage-like structures.⁷³ These voids are filled by a nearby surface atom when the geometry is reoptimised for the Gupta potential, as indicated by the two red arrows in Fig. 5c, thus recovering the ninth isomer from Table 1 (with additional and less significant “breathing” of other atoms). Note that, although the two levels of theory do not exactly agree on the GM structure, they both predict the three high-symmetry structures to be very unfavourable, particularly the cuboctahedron. It is also noteworthy that the DFT energy differences (ΔE_DFT values in Table 1) are an order of magnitude larger than the Gupta energy differences (ΔE_G values), suggesting that the frustrated nature of the Gupta energy landscape may not be preserved at higher levels of theory.

Fig. 5 Au₅₅ GM structure at DFT level from two (a, b) different viewpoints along the (vertically oriented) symmetry plain, with the ball-and-stick representation colour-coded as in Fig. 2. In (c), the DFT structure (black atoms) is superimposed over the ninth isomer (white atoms) from Table 1, with red arrows highlighting the main difference between the two geometries.

To conclude our discussion of the naked “Schmid Au₅₅”, we calculate the HOMO–LUMO energy gap and the partial density of states for isomers 1, 6, 8, 9, 14, Mac, Ino, and Cub in Table 1. The HOMO–LUMO gap (calculated using the ΔSCF method,⁷⁴ comparing the ground state of charge +1 and −1 with the neutral system, including the Makov–Payne correction⁷⁵) ranges from 2.18 to 2.26 eV among these isomers. The total density of states for the three symmetric isomers (particularly Mac) shows pronounced peaks around specific values, while for the more ambiguous structures the density is more broadly distributed. Unfortunately, partial density of states does not immediately reveal any links between the electronic structure and the cluster geometry, but we hope to explore this issue in more detail as a separate study.

B. Size dependence

We now consider other Gupta Au_N clusters in the size range 30 ≤ N ≤ 147, starting with a direct comparison between Au₅₅, Au₈₅, and Au₁₄₇. The corresponding disconnectivity graphs are shown in Fig. 2 and 3, with the heat capacities and motif-decomposed occupation probabilities plotted in Fig. 4.

As for Au₅₅, the ideal Au₁₄₇ Mackay icosahedron³² is energetically unfavourable, with the disconnectivity graph based on 10⁴ lowest-lying minima featuring hardly any traces of the ICO motif. The fraction of AMB minima is also reduced (compared to Au₅₅) and the existence of two structurally homogeneous funnels is apparent: one is predominantly DEC, and the other is TWI. Note that the GM of Au₁₄₇ is a 146-atom Marks decahedron,⁷² with six atoms along the fivefold disclination (i.e. the decahedral spine), and an extra adatom on one of the peripheral (100) facets. The DEC motif is the most populated for k_BT ≤ 40 meV (Fig. 4) and is gradually supplanted by the ICO motif in the range 50 ± 10 meV, over which P^AMB_occ and P^TWI_occ also rise up to 0.09 and 0.04, respectively. The onset of the ICO motif is more abrupt than for Au₅₅, still coinciding with a C_V peak, but the corresponding temperature range is noticeably above the range over which δ reaches the value of 0.3. This apparent mismatch between molecular dynamics and the harmonic superposition approximation may be due to our database of minima under-representing the melted region, but it could also be due to an harmonic effects.

Au₈₅ is in some ways intermediate between Au₅₅ and Au₁₄₇. Its disconnectivity graph (see Fig. 3) shows fairly pronounced FCC, TWI, and DEC funnels, and each funnel exhibits a considerable number of AMB minima. In this particular case the distinction between DEC and AMB minima is marginal, because many AMB structures still exhibit a well-defined decahedral spine, albeit with a higher degree of amorphisation and/or locally (poly)icosahedral order at the surface. As a consequence of our motif definitions, the low-temperature C_V peak for Au₈₅ in Fig. 4 straddles two crossover temperatures. At k_BT = 13 meV the TWI motif is supplanted by the DEC motif as the most populated, and at k_BT = 17 meV it is the AMB motif that starts to dominate. The Lindemann index rises at a noticeably lower temperature than the second C_V peak, similar to Au₁₄₇, showing that the discrepancy is not specific to a particular cluster size. Interestingly, the increase in δ for Au₈₅ and Au₁₄₇ seems to better align with the onset of ICO isomers, when P^ICO_occ ceases to be negligible but does not yet dominate, suggesting that the onset of multiple fivefold disclinations and icosahedral cores can be take as an indicator of melting.

From Fig. 4 it is apparent that the maximal value of P^AMB_occ diminishes with cluster size, which is consistent with Bao et al.²⁶ finding amorphous GM only for N < 55. However, our thermodynamic analysis shows that the AMB motif can also dominate in larger clusters at finite temperatures. To explore this avenue further we systematically analysed Au_N clusters with N = 30–147, accumulating a database of about 10⁴–10⁵ low-lying minima for each N. We also determined two crossover temperatures, T_A and T_I, respectively marking when the AMB and ICO motifs become the most populated. The results are summarised in Table 2, together with an estimate of the melting temperature (T_m) range obtained from molecular dynamics simulations using the Lindemann⁶⁶ index defined in eqn (10).

Table 2 The GM energy, motif, and point group (PG) for Gupta Au_N clusters; two crossover temperatures T_A and T_I, respectively marking when the AMB and ICO motifs become the most probable; and the melting temperature T_m, estimated from MD simulations for selected sizes using the Lindemann criterion. Note that T_A = 0 K when the GM motif is AMB, and T_A and T_I are N/A when the corresponding motif is never the most probable for k_BT ≤ 0.1 eV

N	E GM (eV)	Motif	PG	T A (K)	T I (K)	T m (K)	N	E GM (eV)	Motif	PG	T A (K)	T I (K)	T m (K)
a Low-lying isomers of borderline ICO/AMB motif produce reentrant behaviour and a particularly broad heat capacity peak. b The GM motif is first supplanted by the DEC motif, and then by the AMB motif at temperature TA. c The lowest-lying ICO isomer is borderline and has been reclassified as AMB to eliminate artificial reentrant behaviour.
30	−104.743100	AMB	C 3v	0	460	290 ± 50	60	−213.536346	TWI	C 1	142	446	370 ± 30
31	−108.183974	AMB	C 2	0	402		61	−217.254337	TWI	C 3v	223	585
32	−111.794339	AMB	C 3	0	306		62	−220.855915	TWI	C s	209	504
33^a	−115.404251	ICO	C 1	193	339		63^b	−224.532317	TWI	C s	223	316
34	−119.082484	AMB	C 3	0	295		64	−228.254019	DEC	C 2v	299	397
35	−122.678712	AMB	C s	0	330	260 ± 50	65	−231.868903	DEC	C 2v	220	457	370 ± 50
36	−126.240885	AMB	C 2	0	100	280 ± 50	66	−235.547173	DEC	C s	70	476
37	−129.991492	AMB	C 2v	0	248	280 ± 50	67	−239.158600	DEC	C s	207	555
38	−133.584814	AMB	C s	0	309	300 ± 50	68	−242.838649	DEC	C 2v	153	935
39	−137.184370	FCC	C 4v	12	137	300 ± 50	69	−246.450465	DEC	C 1	172	982
40	−140.789863	AMB	C s	0	35	330 ± 50	70	−250.154102	DEC	C s	100	N/A	370 ± 30
41	−144.403059	AMB	C s	0	262		71	−253.959033	DEC	C 2v	339	N/A
42	−148.023359	AMB	C 1	0	367		72	−257.571776	DEC	C s	239	N/A
43	−151.721379	DEC	C 2v	9	339	320 ± 50	73	−261.253928	DEC	C s	95	N/A
44	−155.322847	AMB	C s	0	260		74	−264.922500	DEC	C 5v	297	N/A
45	−158.916745	DEC	C 2v	9	422	330 ± 20	75	−268.761948	DEC	D 5h	397	N/A	380 ± 30
46	−162.598451	AMB	C 3	0	487		76	−272.372659	DEC	C 2v	337	N/A
47	−166.245363	DEC	C 2v	128	441	340 ± 50	77	−275.982269	DEC	C 2v	123	N/A
48	−169.873475	AMB	C 1	0	471		78	−279.588677	DEC	C 2v	158	N/A
49	−173.562095	DEC	D 5h	111	513		79	−283.417486	FCC	O h	385	N/A	420 ± 30
50	−177.090944	TWI	D 3h	35	520	330 ± 20	80	−287.022207	FCC	C 4v	265	N/A	410 ± 30
51	−180.698557	AMB	C 1	0	457		81	−290.660813	DEC	C 2v	169	N/A
52	−184.431342	AMB	C 2v	0	545		82	−294.269106	DEC	C 2v	232	N/A
53	−188.009110	AMB	C 3v	0	552		83	−297.933795	TWI	C s	276	N/A
54	−191.686569	FCC	C 2v	14	404		84	−301.598992	DEC	C s	181	N/A
55	−195.284851	FCC	C 1	70	501	350 ± 20	85	−305.259030	TWI	C 2v	197	N/A	420 ± 30
56	−199.015100	FCC	D 2h	165	344		86	−308.936312	DEC	C 2v	318	N/A
57^c	−202.624795	TWI	C 2v	176	534	360 ± 50	87	−312.550788	DEC	C 2v	367	N/A
58	−206.216538	TWI	C 2	42	483		88	−316.280622	FCC	C s	327	N/A
59	−209.852603	TWI	C 2v	153	295		89	−319.884977	FCC	C s	443	N/A
120	−434.088354	DEC	C s	624	N/A	470 ± 20	90	−323.613329	FCC	C s	436	N/A	420 ± 20
140	−508.052209	DEC	C s	778	N/A	500 ± 30	95	−341.874115	TWI	C s	309	N/A	420 ± 30
147	−533.942249	DEC	C s	N/A	590	490 ± 30	100	−360.342412	FCC	C 1	343	N/A	450 ± 30

---

Table 2 shows that the GM structures are predominantly AMB for N < 54, then primarily FCC or TWI in the size range 54 ≤ N < 64, and mostly DEC in the range 64 ≤ N < 88, which includes the closed-shell Marks decahedron⁷¹ for N = 75. Note that AMB and DEC motifs are in close competition for N between 42 and 49, as discussed in more detail in section IIIC, while the size range 88 ≤ N ≤ 147 shows re-entrance of the FCC/TWI and DEC motifs. As cluster size increases, the GM motif changes in the sequence AMB → FCC/TWI → DEC → FCC/TWI → DEC → …, which differs from the generally expected^41,76 sequence ICO → DEC → FCC/TWI. Furthermore, although the AMB motif is the GM only for N < 54, we find that T_A < T_m and T_A < T_I for most of the cluster sizes considered, and so amorphous structures still dominate for N ≥ 54 at temperatures between T_A and min{T_m, T_I}. Also note that T_I > T_m or T_I is N/A in many cases, which implies general thermodynamic instability of well-defined icosahedral order in the solid state.

In passing we confirm that the GM of Au₃₈ is a low-symmetry AMB isomer, as discovered by Garzón et al.,²² with the symmetric truncated octahedron (of FCC motif) higher in energy by 6 meV. The ordering is reversed when one extra atom is added: the GM of Au₃₉ is FCC, comprising the 38-atom truncated octahedron with the extra atom placed on one of the (100) facets, essentially converting the facet into a vertex. This behaviour is again consistent with energetically unfavourable (100) facets,³⁵ making Au₃₉ the smallest cluster with a single-crystal fcc GM, which beats the next-lowest AMB isomer by less than 2 meV. Given these small energy differences it is reasonable to expect the ordering to vary for different models and levels of theory.

C. Geometric odd–even behaviour

Interestingly, near degeneracy of the AMB and DEC motifs in the size-range 42 ≤ N ≤ 49 causes the GM structure to alternate, with the DEC motif prevailing for odd N. This behaviour is illustrated in Fig. 6, providing an interesting example of geometric odd–even effect in atomic clusters. The lowest-lying DEC isomer is not the GM for even N, when there is an unpaired adatom on one of the peripheral (100) facets of the incomplete Marks decahedron. Although this odd–even effect is unlikely to be general, and indeed it is not as pronounced in Sutton–Chen clusters⁷⁷ (breaking down for N = 45), it nonetheless illustrates a more extreme case of the reentrant phase behaviour reported for silver clusters.⁷⁸ We also note that in several cases (e.g. N = 43, 47, and 49) the structural differences between the lowest-lying DEC and AMB minima seem fairly minor: the AMB isomer still exhibits a discernible decahedral spine, which is slightly crooked due to more localised fivefold disclinations near the surface.

Fig. 6 The GM (bottom row) and the lowest-lying minimum of a different motif (top row) illustrated for Au_N clusters in the size range exhibiting odd–even behaviour, with the GM alternating between AMB and DEC motifs. Red spheres represent ico atoms and the arrows keep track of the DEC motif with four atoms along the decahedral spine.

To check that the range of odd–even behaviour in Table 2 is not an artefact of our motif definitions, but actually has thermodynamic implications, we consider the heat capacity in the size range 42 ≤ N ≤ 49. Fig. 7 shows two well-defined C_V peaks featuring for odd N, one near the temperatures T_A and the other closer to T_I. The low-temperature peak in each case marks the crossover from DEC to AMB motif, while the high temperature peak straddles the melting range and a crossover from the AMB to the ICO motif as the most populated. The low-temperature peaks are absent for even N, consistent with T_A = 0. The results for odd N = 43, 47, and 49, on the other hand, show that thermally activated manifestation of seemingly minor local icosahedral order near the surface of a cluster can produce a well-defined peak in the heat capacity.

Fig. 7 Heat capacity (C_V) versus temperature (k_BT) for Au_N clusters with even (top) and odd (bottom) N in the range 42 ≤ N ≤ 49. For Au₄₈, the indicated C_V shoulder corresponds to the isomer illustrated in the inset supplanting the GM shown in Fig. 6.

D. Vibrational entropy and fivefold disclinations

Although the C_V plots in Fig. 7 do not exhibit a pronounced low-temperature peak for even N, there are still some identifiable premelting features in the form of a low-temperature shoulder to the main (melting) peak. These features arise from inhomogeneities within the AMB motif. For instance, the low-temperature shoulder for Au₄₈ marks a transition from the GM (see Fig. 6) to a higher energy minimum (see Fig. 7 inset) of the same AMB motif, but with considerably more ico atoms. The highest ico coordination of ico atoms is below seven in both minima, which is why they are both classified as AMB. Hence, there is scope for splitting the AMB motif into finer submotifs, defined by additional criteria on the number and topology of ico atoms, or using a more general approach⁷⁹ based on the temperature derivative of occupation probabilities for individual minima, which can help better interpret various C_V features. In the present model we generally observe that vibrational entropy tends to be larger for minima with more ico atoms, indicative of normal modes softening with increasing fivefold disclination density. This observation explains why thermally activated morphological transitions generally follow the order FCC/TWI → DEC → AMB → ICO, which corresponds to increasing number of ico atoms and increasing vibrational entropy. Note that equilibrium transitions of type FCC/TWI → DEC occur for N = 63, 79, 83, 85, 89 and 90, as indicated by the superscript b in Table 2.

E. Rearrangement mechanisms

We now consider the rearrangement mechanisms for selected morphological transitions, focusing on the fastest transition pathways in Au₅₅, Au₈₅ and Au₁₄₇. We start with the pathway connecting the symmetric cuboctahedral isomer of Au₅₅ with the fcc GM. Given that both endpoints are single-crystal fcc, our intuition suggests a transition pathway based on a sequence of atom-hops on the cluster surface. However, it turns out that the fastest route in and out of the cuboctahedral local minimum is via a particular transition state of point group T_d, corresponding to the diamond–square–diamond mechanism,⁴³ which connects the cuboctahedron to the closed-shell Mackay icosahedron.²⁰ The complete pathway is illustrated in Fig. 8, showing how the intervening Mackay icosahedron rearranges into the fcc GM via a sequence of local minima with no discernible order.

Fig. 8 Potential energy profiles of the fastest discrete pathway between selected endpoints (isomers Cub, 1, 5 and 8 from Table 1). The horizontal and the vertical axes correspond to the (Gupta) potential energy and the discrete path length, respectively, with the same scaling for all the profiles. Solid lines trace pathways for k_BT = 26 meV (room temperature) and dashed for 5 meV. Filled and unfilled symbols correspond to minima and saddles, respectively. The structure of three minima (labelled A, B and C) is shown with the amb and ico atoms highlighted and colour-coded as in Fig. 2.

Although the lowest-lying AMB isomer is only 36 meV higher in energy than the GM, the lowest overall energy barrier is about 0.5 eV, and the corresponding pathway involves minima with partially disordered geometries, such as snapshot B in Fig. 8. Hence, if this transition is observed, it would most likely resemble partial melting followed by crystallisation into the fcc structure. Interestingly, the maximal degree of amorphisation along the fastest path decreases as the temperature increases, which can be seen by comparing the fastest pathways at k_BT = 26 meV and 5 meV. This somewhat counter-intuitive trend is not unexpected considering that most low-lying minima of Au₅₅ exhibit fairly amorphous or ambiguous structure (recall Fig. 2). It is also noteworthy that the number of steps in the fastest pathway is more than halved when the temperature is increased from k_BT = 5 meV to room temperature (26 meV), showing that the fastest mechanism is temperature dependent.

Fig. 8 also shows the potential energy profile for the optimal pathway between isomers five and eight in Table 1 at two different temperatures. The profiles involve considerably lower energy barriers compared to the other pathways, consistent with the disconnectivity graph in Fig. 2, and showing no significant change with temperature.

Recall that the most competitive motifs for N > 56 at low temperatures are DEC and TWI, so we examine if the fastest pathway between these two motifs follows any particular mechanism, using Au₈₅ and Au₁₄₇ at room temperature as our test cases. For Au₈₅, Fig. 9 shows that part of the twin boundary in the lowest-lying TWI isomer is preserved along the pathway to the lowest-lying DEC isomer. However, formation of the fivefold twin axis is accompanied by a significant level of disorder, which would probably be interpreted as partial melting if observed in experiments. A more significant level of disorder occurs in the Au₁₄₇ pathway (see Fig. 9), where (unlike Au₈₅) the amorphous intermediates exhibit multiple well-defined fivefold disclinations, somewhat resembling a distorted icosahedron. Most of these amorphous intermediates are about 1 eV higher in energy than the TWI and DEC endpoints, while among themselves they are separated by relatively small energy barriers, which also suggests that the calculated pathway passes through a liquid-like state.

Fig. 9 Potential energy profiles for the fastest discrete pathway from the lowest-lying TWI minimum (leftmost) to the lowest-lying DEC minimum (rightmost) for Au₈₅ (top) and Au₁₄₇ (bottom) at room temperature (k_BT = 26 meV). The horizontal and the vertical axes, the meaning of filled/unfilled symbols, and the colour-coding of the atoms in the selected snapshots are as in Fig. 8.

IV. Summary and conclusions

We have applied the energy landscapes framework to map the equilibrium morphologies of Au_N clusters in the size range 30 ≤ N ≤ 147, modelled using the Gupta potential, and we examined the rearrangement mechanisms between some competing structures. Our results complement previous theoretical and experimental studies and shed new light on finite-size and temperature effects.

The global minimum (GM) structure of most clusters in the size range 30 ≤ N ≤ 53 was found to be fairly ambiguous, but with discernible fivefold disclinations. We identified a case of geometric odd–even behaviour in the range 42 ≤ N ≤ 49, where the GM alternates between a structure with just one disclination and a structure with multiple disclinations. Global minima of larger clusters (54 ≤ N ≤ 147) are typically a lamellar-twinned or single-crystal fcc lump, or a decahdral motif with a single fivefold disclination. While this structural change is consistent with previous global optimisation studies,²⁶ our thermodynamic analysis further shows that ambiguous motifs with multiple fivefold disclinations can still dominate for N ≤ 85 (except N = 71, 75, 76, and 79) at around room temperature.

In most cases where the GM has at most one disclination, we found thermally-driven morphological transitions below the size-dependent melting temperature. These finite-system analogues of a solid–solid phase transition often coincide with a premelting feature in the heat capacity curve, and they typically correspond to the GM occupation probability dropping below that of multiple ambiguous structures above a system-specific temperature. In certain cases (N = 63, 79, 83, 85, 89, 90) we found two consecutive premelting transitions between distinctly different motifs, with the corresponding crossover temperatures correlating with a smeared shoulder or peak in the heat capacity. In all cases, the higher-energy phase exhibits more fivefold disclinations and higher vibrational entropy.

We calculated the fastest solid–solid transition pathways for Au_N with N = 55, 85, 147, where the energy barrier separating some of the competing motifs is up to 1 eV. While such energy barriers are difficult to treat using conventional simulation methods, which makes direct simulation of the rearrangement mechanism unfeasible, discrete path sampling and harmonic transition state theory provide a useful approximation that performs best at low temperatures. The calculated pathways pass through many metastable intermediates with fivefold disclinations and/or a high degree of amorphisation, consistent with the melt-freeze scenario described by Koga et al.⁷ for larger gold clusters.

Finally, we confirmed a previously reported^26,29,80 fcc GM structure of point group C₁ for Au₅₅, but our thermodynamic analysis revealed that several distorted icosahedra collectively become more favourable at temperatures above 50 Kelvin. Interestingly, the room-temperature occupation probability of a particular isomer was found to be consistent with electron microscopy observations of Wang and Palmer.⁸ This apparent agreement is surprising, because the Gupta potential energy landscape for Au₅₅ clearly disagrees with density functional theory (DFT): relaxing the geometry of sixteen lowest-lying Gupta minima using DFT produced a markedly different energetic ordering, with a new putative GM of point group C_s at the DFT level. While the consistency between our Gupta-level calculations and previous experiments may well be fortuitous, it nonetheless highlights the importance of accounting for thermal fluctuations in geometry—something that is often overlooked when comparing empirical potentials with DFT. In future studies it would be interesting to investigate the equilibrium thermodynamics of Au₅₅ (and other clusters) at the DFT level, where it is also possible to compute the normal-mode frequencies for the harmonic superposition approximation.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was financially supported by EPSRC grant EP/J010847/1 and the ERC. FB is also grateful to the UK Materials and Molecular Modelling Hub for computational resources, which were partially funded by EPSRC grants EP/P020194/1 and EP/J010812/1.

References

1. P. Buffat and J.-P. Borel, Phys. Rev. A, 1976, 13, 2287.
2. A. Sanchez, S. Abbet, U. Heiz, W.-D. Schneider, H. Häkkinen, R. N. Barnett and U. Landman, J. Phys. Chem. A, 1999, 103, 9573; W. Rechberger, A. Hohenau, A. Leitner, J. Krenn, B. Lamprecht and F. Aussenegg, Opt. Commun., 2003, 220, 137; P. Pyykkö, Angew. Chem., Int. Ed., 2004, 43, 4412; M. Walter, J. Akola, O. Lopez-Acevedo, P. D. Jadzinsky, G. Calero, C. J. Ackerson, R. L. Whetten, H. Grönbeck and H. Häkkinen, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 9157.
3. M.-C. Daniel and D. Astruc, Chem. Rev., 2004, 104, 293; S. Eustis and M. A. El-Sayed, Chem. Soc. Rev., 2006, 35, 209; A. S. K. Hashmi and G. J. Hutchings, Angew. Chem., Int. Ed., 2006, 45, 7896; C. Louis and O. Pluchery, Gold nanoparticles for physics, chemistry and biology, World Scientific, 2012.
4. D. Cunningham, W. Vogel, H. Kageyama, S. Tsubota and M. Haruta, J. Catal., 1998, 177, 1.
5. J.-O. Bovin, R. Wallenberg and D. J. Smith, Nature, 1985, 317, 47; D. J. Smith, A. K. Petford-Long, L. R. Wallenberg and J. O. Bovin, Science, 1986, 233, 872; S. Iijima and T. Ichihashi, Phys. Rev. Lett., 1986, 56, 616; M. Mitome, Y. Tanishiro and K. Takayanagi, Z. Phys. D: At., Mol. Clusters, 1989, 12, 45.
6. C. L. Cleveland, U. Landman, T. G. Schaaff, M. N. Shafigullin, P. W. Stephens and R. L. Whetten, Phys. Rev. Lett., 1997, 79, 1873; C. L. Cleveland, U. Landman, M. N. Shafigullin, P. W. Stephens and R. L. Whetten, Z. Phys. D: At., Mol. Clusters, 1997, 40, 503.
7. K. Koga, T. Ikeshoji and K.-i. Sugawara, Phys. Rev. Lett., 2004, 92, 115507.
8. Z. Wang and R. E. Palmer, Nano Lett., 2012, 12, 5510.
9. L. D. Marks, P. Ajayan and J. Dundurs, Ultramicroscopy, 1986, 20, 77; J. Dundurs, L. D. Marks and P. M. Ajayan, Philos. Mag. A, 1988, 57, 605; P. M. Ajayan and L. D. Marks, Phys. Rev. Lett., 1988, 60, 585; P. M. Ajayan and L. D. Marks, Phys. Rev. Lett., 1989, 63, 279; P. M. Ajayan and L. D. Marks, Phase Transitions, 1990, 24–26, 229.
10. W. Krakow, M. Jóse-Yacamán and J. L. Aragón, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 10591.
11. D. J. Wales, Energy Landscapes, Cambridge University Press, Cambridge, 2003.
12. G. Schmid, R. Pfeil, R. Boese, F. Bandermann, S. Meyer, G. H. M. Calis and J. W. A. van der Velden, Chem. Ber., 1981, 114, 3634.
13. G. Schmid, Chem. Soc. Rev., 2008, 37, 1909.
14. M. A. Marcus, M. P. Andrews, J. Zegenhagen, A. S. Bommannavar and P. Montano, Phys. Rev. B: Condens. Matter Mater. Phys., 1990, 42, 3312.
15. W. Vogel, B. Rosner and B. Tesche, J. Phys. Chem., 1993, 97, 11611; D. H. Rapoport, W. Vogel, H. Cölfen and R. Schlögl, J. Phys. Chem. B, 1997, 101, 4175.
16. Y. Pei, N. Shao, Y. Gao and X. C. Zeng, ACS Nano, 2010, 4, 2009.
17. H.-G. Boyen, G. Kästle, F. Weigl, B. Koslowski, C. Dietrich, P. Ziemann, J. P. Spatz, S. Riethmüller, C. Hartmann, M. Möller, G. Schmid, M. G. Garnier and P. Oelhafen, Science, 2002, 297, 1533.
18. M. Turner, V. B. Golovko, O. P. Vaughan, P. Abdulkin, A. Berenguer-Murcia, M. S. Tikhov, B. F. Johnson and R. M. Lambert, Nature, 2008, 454, 981.
19. S. Sawada and S. Sugano, Z. Phys. D: At., Mol. Clusters, 1992, 24, 377.
20. D. J. Wales and L. J. Munro, J. Phys. Chem., 1996, 100, 2053.
21. I. L. Garzón and A. Posada-Amarillas, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11796.
22. I. L. Garzón, K. Michaelian, M. R. Beltrán, A. Posada-Amarillas, P. Ordejón, E. Artacho, D. Sánchez-Portal and J. M. Soler, Phys. Rev. Lett., 1998, 81, 1600.
23. I. L. Garzón, J. A. Reyes-Nava, J. I. Rodríguez-Hernández, I. Sigal, M. R. Beltrán and K. Michaelian, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 66, 073403.
24. X. López-Lozano, L. A. Pérez and I. L. Garzón, Phys. Rev. Lett., 2006, 97, 233401.
25. S. Tehuacanero, R. Herrera, M. Avalos and M. Yacamàn, Acta Metall. Mater., 1992, 40, 1663.
26. K. Bao, S. Goedecker, K. Koga, F. Lançon and A. Neelov, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 041405.
27. J. P. K. Doye and D. J. Wales, New J. Chem., 1998, 22, 733.
28. M. Eckhoff, D. Schebarchov and D. J. Wales, J. Phys. Chem. Lett., 2017, 8, 5402.
29. B. E. Husic, D. Schebarchov and D. J. Wales, Nanoscale, 2016, 8, 18326.
30. F. Cleri and V. Rosato, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 48, 22.
31. H. Häkkinen, M. Moseler, O. Kostko, N. Morgner, M. A. Hoffmann and B. v. Issendorff, Phys. Rev. Lett., 2004, 93, 093401.
32. A. L. Mackay, Acta Crystallogr., 1962, 15, 916.
33. S. Ino, J. Phys. Soc. Jpn., 1969, 27, 941.
34. J. P. K. Doye, D. J. Wales and R. S. Berry, J. Chem. Phys., 1995, 103, 4234; J. M. Soler, M. R. Beltrán, K. Michaelian, I. L. Garzón, P. Ordejón, D. Sánchez-Portal and E. Artacho, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 61, 5771.
35. F. Baletto, R. Ferrando, A. Fortunelli, F. Montalenti and C. Mottet, J. Chem. Phys., 2002, 116, 3856; W. Huang, M. Ji, C.-D. Dong, X. Gu, L.-M. Wang, X. G. Gong and L.-S. Wang, ACS Nano, 2008, 2, 897.
36. N. Tarrat, M. Rapacioli, J. Cuny, J. Morillo, J.-L. Heully and F. Spiegelman, Comput. Theor. Chem., 2017, 1107, 102 , structure prediction of nanoclusters from global optimization techniques: computational strategies. X. Wu and Y. Sun, J. Nanopart. Res., 2017, 19, 201.
37. C. L. Cleveland, W. D. Luedtke and U. Landman, Phys. Rev. Lett., 1998, 81, 2036; C. L. Cleveland, W. D. Luedtke and U. Landman, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 60, 5065.
38. D. Schebarchov and S. C. Hendy, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 121402; D. Schebarchov and S. C. Hendy, Eur. Phys. J. D, 2007, 43, 11.
39. G. A. Breaux, C. M. Neal, B. Cao and M. F. Jarrold, Phys. Rev. Lett., 2005, 94, 173401; K. G. Steenbergen and N. Gaston, Chem. – Eur. J., 2015, 21, 2862.
40. V. A. Mandelshtam and P. A. Frantsuzov, J. Chem. Phys., 2006, 124, 204511; V. A. Sharapov, D. Meluzzi and V. A. Mandelshtam, Phys. Rev. Lett., 2007, 98, 105701; V. A. Sharapov and V. A. Mandelshtam, J. Phys. Chem. A, 2007, 111, 10284.
41. J. P. K. Doye and F. Calvo, Phys. Rev. Lett., 2001, 86, 3570; J. P. K. Doye and F. Calvo, J. Chem. Phys., 2002, 116, 8307.
42. F. H. Stillinger and T. A. Weber, Phys. Rev. A, 1982, 25, 978; D. J. Wales, Mol. Phys., 1993, 78, 151; G. Franke, E. R. Hilf and P. Borrmann, J. Chem. Phys., 1993, 98, 3496; F. Calvo, J. P. K. Doye and D. J. Wales, J. Chem. Phys., 2001, 115, 9627.
43. J. Uppenbrink and D. J. Wales, J. Chem. Soc., Faraday Trans., 1991, 87, 215; D. J. Wales and J. Uppenbrink, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 12342; A. A. Tal, E. P. Munger and I. A. Abrikosov, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 020102.
44. K. Rossi and F. Baletto, Phys. Chem. Chem. Phys., 2017, 19, 11057; K. Rossi, Y. Soon, L. Pavan and F. Baletto, Eur. Phys. J. B, 2017 , arXiv:1702.07088.
45. R. P. Gupta, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 23, 6265.
46. V. Rosato, M. Guillope and B. Legrand, Philos. Mag. A, 1989, 59, 321.
47. J.-P. Palomares-Baez, E. Panizon and R. Ferrando, Nano Lett., 2017, 17, 5394.
48. Z. Li and H. A. Scheraga, Proc. Natl. Acad. Sci. U. S. A., 1987, 84, 6611.
49. D. J. Wales and J. P. K. Doye, J. Phys. Chem. A, 1997, 101, 5111.
50. D. Schebarchov and D. J. Wales, Phys. Chem. Chem. Phys., 2015, 17, 28331.
51. D. J. Wales et al. , GMIN: A program for basin-hopping global optimisation, basin-sampling, and parallel tempering, http://www-wales.ch.cam.ac.uk/software.html.
52. J. M. Carr, S. A. Trygubenko and D. J. Wales, J. Chem. Phys., 2005, 122, 234903.
53. D. J. Wales et al. , OPTIM: A program for geometry optimisation and pathway calculations, http://www-wales.ch.cam.ac.uk/software.html.
54. S. A. Trygubenko and D. J. Wales, J. Chem. Phys., 2004, 120, 2082.
55. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901; G. Henkelman and H. Jónsson, J. Chem. Phys., 2000, 113, 9978.
56. L. J. Munro and D. J. Wales, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 3969.
57. J. Nocedal, Math. Comput., 1980, 35, 773; D. Liu and J. Nocedal, Math. Prog., 1989, 45, 503.
58. D. J. Wales, Mol. Phys., 2002, 100, 3285.
59. J. M. Carr and D. J. Wales, Phys. Chem. Chem. Phys., 2009, 11, 3341.
60. D. J. Wales et al. , PATHSAMPLE: A program for generating connected stationary point databases and extracting global kinetics, http://www-wales.ch.cam.ac.uk/software.html.
61. O. M. Becker and M. Karplus, J. Chem. Phys., 1997, 106, 1495.
62. G. H. Vineyard, J. Phys. Chem. Solids, 1957, 3, 121.
63. J. D. Honeycutt and H. C. Andersen, J. Phys. Chem., 1987, 91, 4950; D. Faken and H. Jónsson, Comput. Mater. Sci., 1994, 2, 279.
64. S. C. Hendy and J. P. K. Doye, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 66, 235402.
65. M. Paterlini and D. M. Ferguson, Chem. Phys., 1998, 236, 243.
66. F. A. Lindemann, Phys. Z., 1910, 11, 609; Y. Zhou, M. Karplus, K. D. Ball and R. S. Berry, J. Chem. Phys., 2002, 116, 2323.
67. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502.
68. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
69. M. J. Piotrowski, C. G. Ungureanu, P. Tereshchuk, K. E. A. Batista, A. S. Chaves, D. Guedes-Sobrinho and J. L. F. Da Silva, J. Phys. Chem. C, 2016, 120, 28844.
70. E. Aprà, F. Baletto, R. Ferrando and A. Fortunelli, Phys. Rev. Lett., 2004, 93, 065502.
71. L. D. Marks, Philos. Mag. A, 1984, 49, 81.
72. T. Ben-David, Y. Lereah, G. Deutscher, J. M. Penisson, A. Bourret, R. Kofman and P. Cheyssac, Phys. Rev. Lett., 1997, 78, 2585.
73. J. Wang, J. Jellinek, J. Zhao, Z. Chen, R. B. King and P. von Ragu Schleyer, J. Phys. Chem. A, 2005, 109, 9265; S. Bulusu, X. Li, L.-S. Wang and X. C. Zeng, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 8326; L. Trombach, S. Rampino, L.-S. Wang and P. Schwerdtfeger, Chem. – Eur. J., 2016, 22, 8823.
74. J. Xian, S. Baroni and P. Umari, J. Chem. Phys., 2014, 140, 124101; F. Baletto and R. Ferrando, Phys. Chem. Chem. Phys., 2015, 17, 28256.
75. G. Makov and M. C. Payne, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 51, 4014.
76. A. S. Barnard, N. P. Young, A. I. Kirkland, M. A. van Huis and H. Xu, ACS Nano, 2009, 3, 1431; J. M. Rahm and P. Erhart, Nano Lett., 2017, 17, 5775.
77. J. P. K. Doye and D. J. Wales, J. Chem. Phys., 2002, 116, 3777.
78. F. Baletto, C. Mottet and R. Ferrando, Phys. Rev. Lett., 2000, 84, 5544.
79. D. J. Wales, Phys. Rev. E, 2017, 95, 030105(R).
80. J. P. K. Doye and D. J. Wales, Phys. Rev. Lett., 1998, 80, 1357.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c7nr07123j

---