Open Access Article
This Open Access Article is licensed under a
Creative Commons Attribution 3.0 Unported Licence

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

D. Schebarchov *a, F. Baletto b and D. J. Wales a
aUniversity Chemical Laboratories, Lensfield Road, Cambridge CB2 1EW, UK. E-mail: Dmitri.Schebarchov@gmail.com; dw34@cam.ac.uk
bDepartment of Physics, King's College London, London WC2R 2LS, UK. E-mail: francesca.baletto@kcl.ac.uk

Received 24th September 2017 , Accepted 15th December 2017

First published on 15th December 2017


We consider finite-size and temperature effects on the structure of model AuN 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.3 Some novel properties, such as the unexpected catalytic activity,4 have been linked to specific cluster morphologies, and the perceived structure–property relationship has motivated many studies5–8 of the atomic-level structure evolution. However, the diversity of morphologies and quasimelting9,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 framework11 to shed new light on the polymorphism in model gold clusters.

One of the first and still widely studied gold nanoparticles is the “Schmid Au55” cluster,12,13 formulated as Au55 (PPh3)12 Cl6. The Au55 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.,15 who found that a model with icosahedral Au55 core produced a better fit to the available X-ray powder diffraction data. More recently, Pei et al.16 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 Au55 core geometry has still not reached a definitive resolution, exacerbating interpretation of the unusual oxidation resistance17 and high catalytic activity18 of the naked Au55 cluster derived from Au55 (PPh3)12 Cl6. The intriguing chemistry of the naked Au55 may well be at least in part due to a particular geometry,17 but the geometry is unlikely to be cuboctahedral, because previous theoretical calculations19,20 show that the expected lifetime of bare Au55 cuboctahedra is too short to be observed under an electron microscope, even in the presence of a substrate.19


image file: c7nr07123j-f1.tif
Fig. 1 The structure of some Au55 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 Au55, with “amorphous” signifying distorted icosahedral order (see Fig. 1b and c), which can exhibit a degree of chirality23 and enantioselective properties.24 The lowest-lying amorphous isomer (Fig. 1b) was proposed as the global minimum (GM). Indeed, distorted25 or amorphized10 icosahedra and the partially chiral8 Au55 structures have been observed under an electron microscope. However, using the same empirical model for Au55 as Garzón et al.,21 Bao et al.26 found a lower energy structure with fcc order, shown in Fig. 1d, which had previously27 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,28 the electron beam, or some other factors. Recently,29 we performed global optimisation of Au55 for a different parametrisation30 of the model used by Garzón et al.21 and Bao et al.,26 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 Palmer8 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 Au55 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 amorphous22 and fcc31 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 Au55 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 Au55 at DFT level. We also verify that high-symmetry morphologies (such as the Mackay icosahedron,32 Ino decahedron,33 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 potential21,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.36 However, the GM alone is not sufficient to explain the finite-temperature behaviour and morphological changes observed in experiments5,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-dependent1 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 landscapes11 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 AuN clusters (30 ≤ N ≤ 147), focusing on their equilibrium thermodyamics in the harmonic superposition approximation,42 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.7

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 Au55 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 AuN 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 Au55 structures reoptimised using DFT.

II. Model and methods

We model the binding in AuN clusters of size N using the Gupta potential30,45,46
 
image file: c7nr07123j-t1.tif(1)
where rij 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 r0 = 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,47 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-hopping48,49 global optimisation aided by systematic surface refinement,50 with atom-vacancy swap candidates scanned in parallel (as implemented in GMIN51). More than 104 lowest-lying minima accumulated during global optimisation provide the input for a harmonic superposition analysis of the thermodynamics.42 This analysis involves writing the canonical partition function as
 
image file: c7nr07123j-t2.tif(2)
where the sum runs over the local minima, Em is the potential energy of minimum m, and kB is the Boltzmann constant. The degeneracy factor11,29
 
image file: c7nr07123j-t3.tif(3)
subsumes entropic contributions, represented by the point group order (om) 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, and
 
image file: c7nr07123j-t4.tif(4)
where λ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
 
pm(T) = gmeEm/kBT/Z(T),(5)
and the vibrational contribution to the heat capacity is
 
image file: c7nr07123j-t5.tif(6)
where image file: c7nr07123j-t6.tif, and the constant kinetic contribution has been omitted. Since the pm are additive, we can define collective occupation probabilities image file: c7nr07123j-t7.tif, where image file: c7nr07123j-t8.tif 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 approach52 implemented in OPTIM.53 The pairwise connection attempts are performed with the doubly-nudged54 elastic band algorithm,55 with the intervening transition state candidates tightly converged using hybrid eigenvector following,56 and intervening local minima converged using the limited-memory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) algorithm of Liu and Nocedal.57 The resulting database of stationary points is further extended using discrete path sampling,58 exploiting previously described59 strategies implemented in PATHSAMPLE.60 We visualise these networks in the form of disconnectivity graphs,61 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.62 That is, the escape rate from a minimum m via an adjacent transition state s is

 
image file: c7nr07123j-t9.tif(7)
where the pre-exponential frequency factor is given by
 
image file: c7nr07123j-t10.tif(8)
with 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 as
 
image file: c7nr07123j-t11.tif(9)
where the sum is over all the directly connected transition states. Note that ν0m→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),63 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 Doye64) 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.65 We used a damping coefficient of 5 ps−1 (equivalent to a time constant of 20 fs), and the Verlet integration time step was set to 10 fs. Each simulation ran for 2 × 107 time steps or more, with the first 4 × 106 steps treated as equilibration and discarded. Over the remaining 160 ns we computed the time- and atom-averaged Lindemann index66
 
image file: c7nr07123j-t12.tif(10)
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 Au55 minima were reoptimised using the Quantum Espresso package67—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.68 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 5d106s1. 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 Au55 minima.69

III. Results and discussion

A. The naked “Schmid Au55

The disconnectivity graph for 500 lowest-lying minima of the Au55 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-like70 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.29
image file: c7nr07123j-f2.tif
Fig. 2 Colour-coded disconnectivity graph for the 500 lowest-lying minima of the Au55 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.,26 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.8 The lowest-lying DEC minimum is fourteenth overall, with the fivefold disclination significantly off centre—in contrast to the Ino33 and Marks71 decahedra, but reminiscent of some pentagonally twinned structures reported for lead nanoparticles.72

Table 1 Potential energy relative to the GM (i.e. ΔE = EEGM), motif, point group (PG), the geometric mean normal mode vibrational frequencies ([v with combining macron]) in terahertz, and the number of distinct adjacent transition states (Nt.s.) for the sixteen lowest-lying minima and three high-symmetry structures of Gupta Au55. Each minimum's occupation probability (pm) and lifetime (τm) in seconds is calculated at room temperature (kBT = 26 meV) using a database of more than 3 × 105 minima and 4 × 105 transition states. The relative energy (ΔEDFT) of each isomer reoptimised at the DFT level is also given, and the chiral isomer (ranked eighth) imaged by Wang and Palmer8 is marked by a dagger
  ΔEG (eV) Motif PG [v with combining macron] (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 Au55 is fcc with point group C1, 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 known35 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 previous21,37 reports of distorted icosahedral order in Au55, snapshots of the AMB minima in Fig. 2 exhibit triangular close-packed facets and multiple fivefold disclinations, which outline the tetrahedral units expected in Mackay32 icosahedra. The ideal Mackay icosahedron, on the other hand, is not even in the lowest-lying 105 minima (in a database of more around 3 × 105 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 CV and the collective occupation probability PXXXocc of each motif (XXX) as a function of kBT > 0, plotted in Fig. 4.


image file: c7nr07123j-f3.tif
Fig. 3 Disconnectivity graphs for 103 minima of Au85 (top) and 104 minima of Au147 (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.

image file: c7nr07123j-f4.tif
Fig. 4 Motif occupation probabilities Pocc (stacked on top of each other), vibrational heat capacity CV, and time-averaged Lindemann index δ plotted against temperature for Au55 (top), Au85 (middle), and Au147 (bottom). Five motif shares (FCC, AMB, ICO, TWI and DEC) of Pocc are colour-coded as in Fig. 2 and 3. The chiral (eighth in Table 1) Au55 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 CV is represented by a thick black line, while thinner black lines correspond to subsets of 10, 102, 103, 104, and 105 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 Au55, PFCCocc is the highest below 70 K (kBT < 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. PAMBocc rapidly grows for kBT 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 observed8 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.29 In both cases, when PAMBocc supplants PFCCocc 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 Au55 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.8 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 Palmer8 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.9

Fig. 4 shows that the ICO motif hardly features in Au55 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 Au55. The rise of PICOocc with kBT in the range 30 meV < kBT < 80 meV, with PAMBocc 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 Au55 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 Cs 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.73 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 (ΔEDFT values in Table 1) are an order of magnitude larger than the Gupta energy differences (ΔEG values), suggesting that the frustrated nature of the Gupta energy landscape may not be preserved at higher levels of theory.


image file: c7nr07123j-f5.tif
Fig. 5 Au55 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 Au55”, 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,74 comparing the ground state of charge +1 and −1 with the neutral system, including the Makov–Payne correction75) 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 AuN clusters in the size range 30 ≤ N ≤ 147, starting with a direct comparison between Au55, Au85, and Au147. 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 Au55, the ideal Au147 Mackay icosahedron32 is energetically unfavourable, with the disconnectivity graph based on 104 lowest-lying minima featuring hardly any traces of the ICO motif. The fraction of AMB minima is also reduced (compared to Au55) and the existence of two structurally homogeneous funnels is apparent: one is predominantly DEC, and the other is TWI. Note that the GM of Au147 is a 146-atom Marks decahedron,72 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 kBT ≤ 40 meV (Fig. 4) and is gradually supplanted by the ICO motif in the range 50 ± 10 meV, over which PAMBocc and PTWIocc also rise up to 0.09 and 0.04, respectively. The onset of the ICO motif is more abrupt than for Au55, still coinciding with a CV 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.

Au85 is in some ways intermediate between Au55 and Au147. 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 CV peak for Au85 in Fig. 4 straddles two crossover temperatures. At kBT = 13 meV the TWI motif is supplanted by the DEC motif as the most populated, and at kBT = 17 meV it is the AMB motif that starts to dominate. The Lindemann index rises at a noticeably lower temperature than the second CV peak, similar to Au147, showing that the discrepancy is not specific to a particular cluster size. Interestingly, the increase in δ for Au85 and Au147 seems to better align with the onset of ICO isomers, when PICOocc 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 PAMBocc diminishes with cluster size, which is consistent with Bao et al.26 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 AuN clusters with N = 30–147, accumulating a database of about 104–105 low-lying minima for each N. We also determined two crossover temperatures, TA and TI, 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 (Tm) range obtained from molecular dynamics simulations using the Lindemann66 index defined in eqn (10).

Table 2 The GM energy, motif, and point group (PG) for Gupta AuN clusters; two crossover temperatures TA and TI, respectively marking when the AMB and ICO motifs become the most probable; and the melting temperature Tm, estimated from MD simulations for selected sizes using the Lindemann criterion. Note that TA = 0 K when the GM motif is AMB, and TA and TI are N/A when the corresponding motif is never the most probable for kBT ≤ 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
33a −115.404251 ICO C 1 193 339 63b −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
57c −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 decahedron71 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 expected41,76 sequence ICO → DEC → FCC/TWI. Furthermore, although the AMB motif is the GM only for N < 54, we find that TA < Tm and TA < TI for most of the cluster sizes considered, and so amorphous structures still dominate for N ≥ 54 at temperatures between TA and min{Tm, TI}. Also note that TI > Tm or TI 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 Au38 is a low-symmetry AMB isomer, as discovered by Garzón et al.,22 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 Au39 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,35 making Au39 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 clusters77 (breaking down for N = 45), it nonetheless illustrates a more extreme case of the reentrant phase behaviour reported for silver clusters.78 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.
image file: c7nr07123j-f6.tif
Fig. 6 The GM (bottom row) and the lowest-lying minimum of a different motif (top row) illustrated for AuN 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 CV peaks featuring for odd N, one near the temperatures TA and the other closer to TI. 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 TA = 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.


image file: c7nr07123j-f7.tif
Fig. 7 Heat capacity (CV) versus temperature (kBT) for AuN clusters with even (top) and odd (bottom) N in the range 42 ≤ N ≤ 49. For Au48, the indicated CV shoulder corresponds to the isomer illustrated in the inset supplanting the GM shown in Fig. 6.

D. Vibrational entropy and fivefold disclinations

Although the CV 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 Au48 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 approach79 based on the temperature derivative of occupation probabilities for individual minima, which can help better interpret various CV 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 Au55, Au85 and Au147. We start with the pathway connecting the symmetric cuboctahedral isomer of Au55 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 Td, corresponding to the diamond–square–diamond mechanism,43 which connects the cuboctahedron to the closed-shell Mackay icosahedron.20 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.
image file: c7nr07123j-f8.tif
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 kBT = 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 kBT = 26 meV and 5 meV. This somewhat counter-intuitive trend is not unexpected considering that most low-lying minima of Au55 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 kBT = 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 Au85 and Au147 at room temperature as our test cases. For Au85, 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 Au147 pathway (see Fig. 9), where (unlike Au85) 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.


image file: c7nr07123j-f9.tif
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 Au85 (top) and Au147 (bottom) at room temperature (kBT = 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 AuN 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,26 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 AuN 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.7 for larger gold clusters.

Finally, we confirmed a previously reported26,29,80 fcc GM structure of point group C1 for Au55, 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.8 This apparent agreement is surprising, because the Gupta potential energy landscape for Au55 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 Cs 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 Au55 (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 CrossRef CAS .
  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 CrossRef CAS ; W. Rechberger, A. Hohenau, A. Leitner, J. Krenn, B. Lamprecht and F. Aussenegg, Opt. Commun., 2003, 220, 137 CrossRef ; P. Pyykkö, Angew. Chem., Int. Ed., 2004, 43, 4412 CrossRef PubMed ; 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 CrossRef PubMed .
  3. M.-C. Daniel and D. Astruc, Chem. Rev., 2004, 104, 293 CrossRef CAS PubMed ; S. Eustis and M. A. El-Sayed, Chem. Soc. Rev., 2006, 35, 209 RSC ; A. S. K. Hashmi and G. J. Hutchings, Angew. Chem., Int. Ed., 2006, 45, 7896 CrossRef PubMed ; C. Louis and O. Pluchery, Gold nanoparticles for physics, chemistry and biology, World Scientific, 2012 Search PubMed .
  4. D. Cunningham, W. Vogel, H. Kageyama, S. Tsubota and M. Haruta, J. Catal., 1998, 177, 1 CrossRef CAS .
  5. J.-O. Bovin, R. Wallenberg and D. J. Smith, Nature, 1985, 317, 47 CrossRef CAS ; D. J. Smith, A. K. Petford-Long, L. R. Wallenberg and J. O. Bovin, Science, 1986, 233, 872 Search PubMed ; S. Iijima and T. Ichihashi, Phys. Rev. Lett., 1986, 56, 616 CrossRef PubMed ; M. Mitome, Y. Tanishiro and K. Takayanagi, Z. Phys. D: At., Mol. Clusters, 1989, 12, 45 CrossRef .
  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 CrossRef CAS ; C. L. Cleveland, U. Landman, M. N. Shafigullin, P. W. Stephens and R. L. Whetten, Z. Phys. D: At., Mol. Clusters, 1997, 40, 503 CrossRef .
  7. K. Koga, T. Ikeshoji and K.-i. Sugawara, Phys. Rev. Lett., 2004, 92, 115507 CrossRef PubMed .
  8. Z. Wang and R. E. Palmer, Nano Lett., 2012, 12, 5510 CrossRef CAS PubMed .
  9. L. D. Marks, P. Ajayan and J. Dundurs, Ultramicroscopy, 1986, 20, 77 CrossRef CAS ; J. Dundurs, L. D. Marks and P. M. Ajayan, Philos. Mag. A, 1988, 57, 605 CrossRef ; P. M. Ajayan and L. D. Marks, Phys. Rev. Lett., 1988, 60, 585 CrossRef PubMed ; P. M. Ajayan and L. D. Marks, Phys. Rev. Lett., 1989, 63, 279 CrossRef PubMed ; P. M. Ajayan and L. D. Marks, Phase Transitions, 1990, 24–26, 229 CrossRef .
  10. W. Krakow, M. Jóse-Yacamán and J. L. Aragón, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 10591 CrossRef CAS .
  11. D. J. Wales, Energy Landscapes, Cambridge University Press, Cambridge, 2003 Search PubMed .
  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 CrossRef CAS .
  13. G. Schmid, Chem. Soc. Rev., 2008, 37, 1909 RSC .
  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 CrossRef CAS .
  15. W. Vogel, B. Rosner and B. Tesche, J. Phys. Chem., 1993, 97, 11611 CrossRef CAS ; D. H. Rapoport, W. Vogel, H. Cölfen and R. Schlögl, J. Phys. Chem. B, 1997, 101, 4175 CrossRef .
  16. Y. Pei, N. Shao, Y. Gao and X. C. Zeng, ACS Nano, 2010, 4, 2009 CrossRef PubMed .
  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 CrossRef CAS PubMed .
  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 CrossRef CAS PubMed .
  19. S. Sawada and S. Sugano, Z. Phys. D: At., Mol. Clusters, 1992, 24, 377 CrossRef CAS .
  20. D. J. Wales and L. J. Munro, J. Phys. Chem., 1996, 100, 2053 CrossRef CAS .
  21. I. L. Garzón and A. Posada-Amarillas, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11796 CrossRef .
  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 CrossRef .
  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 CrossRef .
  24. X. López-Lozano, L. A. Pérez and I. L. Garzón, Phys. Rev. Lett., 2006, 97, 233401 CrossRef PubMed .
  25. S. Tehuacanero, R. Herrera, M. Avalos and M. Yacamàn, Acta Metall. Mater., 1992, 40, 1663 CrossRef CAS .
  26. K. Bao, S. Goedecker, K. Koga, F. Lançon and A. Neelov, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 041405 CrossRef .
  27. J. P. K. Doye and D. J. Wales, New J. Chem., 1998, 22, 733 RSC .
  28. M. Eckhoff, D. Schebarchov and D. J. Wales, J. Phys. Chem. Lett., 2017, 8, 5402 CrossRef CAS PubMed .
  29. B. E. Husic, D. Schebarchov and D. J. Wales, Nanoscale, 2016, 8, 18326 RSC .
  30. F. Cleri and V. Rosato, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 48, 22 CrossRef CAS .
  31. H. Häkkinen, M. Moseler, O. Kostko, N. Morgner, M. A. Hoffmann and B. v. Issendorff, Phys. Rev. Lett., 2004, 93, 093401 CrossRef PubMed .
  32. A. L. Mackay, Acta Crystallogr., 1962, 15, 916 CrossRef CAS .
  33. S. Ino, J. Phys. Soc. Jpn., 1969, 27, 941 CrossRef CAS .
  34. J. P. K. Doye, D. J. Wales and R. S. Berry, J. Chem. Phys., 1995, 103, 4234 CrossRef CAS ; 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 CrossRef .
  35. F. Baletto, R. Ferrando, A. Fortunelli, F. Montalenti and C. Mottet, J. Chem. Phys., 2002, 116, 3856 CrossRef CAS ; W. Huang, M. Ji, C.-D. Dong, X. Gu, L.-M. Wang, X. G. Gong and L.-S. Wang, ACS Nano, 2008, 2, 897 CrossRef PubMed .
  36. N. Tarrat, M. Rapacioli, J. Cuny, J. Morillo, J.-L. Heully and F. Spiegelman, Comput. Theor. Chem., 2017, 1107, 102 CrossRef CAS  , structure prediction of nanoclusters from global optimization techniques: computational strategies. X. Wu and Y. Sun, J. Nanopart. Res., 2017, 19, 201 CrossRef .
  37. C. L. Cleveland, W. D. Luedtke and U. Landman, Phys. Rev. Lett., 1998, 81, 2036 CrossRef CAS ; C. L. Cleveland, W. D. Luedtke and U. Landman, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 60, 5065 CrossRef .
  38. D. Schebarchov and S. C. Hendy, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 121402 CrossRef CAS ; D. Schebarchov and S. C. Hendy, Eur. Phys. J. D, 2007, 43, 11 CrossRef .
  39. G. A. Breaux, C. M. Neal, B. Cao and M. F. Jarrold, Phys. Rev. Lett., 2005, 94, 173401 CrossRef CAS PubMed ; K. G. Steenbergen and N. Gaston, Chem. – Eur. J., 2015, 21, 2862 CrossRef PubMed .
  40. V. A. Mandelshtam and P. A. Frantsuzov, J. Chem. Phys., 2006, 124, 204511 CrossRef CAS PubMed ; V. A. Sharapov, D. Meluzzi and V. A. Mandelshtam, Phys. Rev. Lett., 2007, 98, 105701 CrossRef PubMed ; V. A. Sharapov and V. A. Mandelshtam, J. Phys. Chem. A, 2007, 111, 10284 CrossRef PubMed .
  41. J. P. K. Doye and F. Calvo, Phys. Rev. Lett., 2001, 86, 3570 CrossRef CAS PubMed ; J. P. K. Doye and F. Calvo, J. Chem. Phys., 2002, 116, 8307 CrossRef .
  42. F. H. Stillinger and T. A. Weber, Phys. Rev. A, 1982, 25, 978 CrossRef CAS ; D. J. Wales, Mol. Phys., 1993, 78, 151 CrossRef ; G. Franke, E. R. Hilf and P. Borrmann, J. Chem. Phys., 1993, 98, 3496 CrossRef ; F. Calvo, J. P. K. Doye and D. J. Wales, J. Chem. Phys., 2001, 115, 9627 CrossRef .
  43. J. Uppenbrink and D. J. Wales, J. Chem. Soc., Faraday Trans., 1991, 87, 215 RSC ; D. J. Wales and J. Uppenbrink, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 12342 CrossRef CAS ; A. A. Tal, E. P. Munger and I. A. Abrikosov, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 020102 CrossRef .
  44. K. Rossi and F. Baletto, Phys. Chem. Chem. Phys., 2017, 19, 11057 RSC ; K. Rossi, Y. Soon, L. Pavan and F. Baletto, Eur. Phys. J. B, 2017 Search PubMed  , arXiv:1702.07088.
  45. R. P. Gupta, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 23, 6265 CrossRef CAS .
  46. V. Rosato, M. Guillope and B. Legrand, Philos. Mag. A, 1989, 59, 321 CrossRef .
  47. J.-P. Palomares-Baez, E. Panizon and R. Ferrando, Nano Lett., 2017, 17, 5394 CrossRef CAS PubMed .
  48. Z. Li and H. A. Scheraga, Proc. Natl. Acad. Sci. U. S. A., 1987, 84, 6611 CrossRef CAS .
  49. D. J. Wales and J. P. K. Doye, J. Phys. Chem. A, 1997, 101, 5111 CrossRef CAS .
  50. D. Schebarchov and D. J. Wales, Phys. Chem. Chem. Phys., 2015, 17, 28331 RSC .
  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 Search PubMed .
  52. J. M. Carr, S. A. Trygubenko and D. J. Wales, J. Chem. Phys., 2005, 122, 234903 CrossRef PubMed .
  53. D. J. Wales et al. , OPTIM: A program for geometry optimisation and pathway calculations, http://www-wales.ch.cam.ac.uk/software.html Search PubMed .
  54. S. A. Trygubenko and D. J. Wales, J. Chem. Phys., 2004, 120, 2082 CrossRef CAS PubMed .
  55. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901 CrossRef CAS ; G. Henkelman and H. Jónsson, J. Chem. Phys., 2000, 113, 9978 CrossRef .
  56. L. J. Munro and D. J. Wales, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 3969 CrossRef CAS .
  57. J. Nocedal, Math. Comput., 1980, 35, 773 CrossRef ; D. Liu and J. Nocedal, Math. Prog., 1989, 45, 503 CrossRef .
  58. D. J. Wales, Mol. Phys., 2002, 100, 3285 CrossRef CAS .
  59. J. M. Carr and D. J. Wales, Phys. Chem. Chem. Phys., 2009, 11, 3341 RSC .
  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 Search PubMed .
  61. O. M. Becker and M. Karplus, J. Chem. Phys., 1997, 106, 1495 CrossRef CAS .
  62. G. H. Vineyard, J. Phys. Chem. Solids, 1957, 3, 121 CrossRef CAS .
  63. J. D. Honeycutt and H. C. Andersen, J. Phys. Chem., 1987, 91, 4950 CrossRef CAS ; D. Faken and H. Jónsson, Comput. Mater. Sci., 1994, 2, 279 CrossRef .
  64. S. C. Hendy and J. P. K. Doye, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 66, 235402 CrossRef .
  65. M. Paterlini and D. M. Ferguson, Chem. Phys., 1998, 236, 243 CrossRef CAS .
  66. F. A. Lindemann, Phys. Z., 1910, 11, 609 Search PubMed ; Y. Zhou, M. Karplus, K. D. Ball and R. S. Berry, J. Chem. Phys., 2002, 116, 2323 CrossRef CAS .
  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 CrossRef PubMed .
  68. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865 CrossRef CAS PubMed .
  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 CAS .
  70. E. Aprà, F. Baletto, R. Ferrando and A. Fortunelli, Phys. Rev. Lett., 2004, 93, 065502 CrossRef PubMed .
  71. L. D. Marks, Philos. Mag. A, 1984, 49, 81 CrossRef CAS .
  72. T. Ben-David, Y. Lereah, G. Deutscher, J. M. Penisson, A. Bourret, R. Kofman and P. Cheyssac, Phys. Rev. Lett., 1997, 78, 2585 CrossRef .
  73. J. Wang, J. Jellinek, J. Zhao, Z. Chen, R. B. King and P. von Ragu Schleyer, J. Phys. Chem. A, 2005, 109, 9265 CrossRef CAS PubMed ; S. Bulusu, X. Li, L.-S. Wang and X. C. Zeng, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 8326 CrossRef PubMed ; L. Trombach, S. Rampino, L.-S. Wang and P. Schwerdtfeger, Chem. – Eur. J., 2016, 22, 8823 CrossRef PubMed .
  74. J. Xian, S. Baroni and P. Umari, J. Chem. Phys., 2014, 140, 124101 CrossRef CAS PubMed ; F. Baletto and R. Ferrando, Phys. Chem. Chem. Phys., 2015, 17, 28256 RSC .
  75. G. Makov and M. C. Payne, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 51, 4014 CrossRef CAS .
  76. A. S. Barnard, N. P. Young, A. I. Kirkland, M. A. van Huis and H. Xu, ACS Nano, 2009, 3, 1431 CrossRef CAS PubMed ; J. M. Rahm and P. Erhart, Nano Lett., 2017, 17, 5775 CrossRef PubMed .
  77. J. P. K. Doye and D. J. Wales, J. Chem. Phys., 2002, 116, 3777 CrossRef CAS .
  78. F. Baletto, C. Mottet and R. Ferrando, Phys. Rev. Lett., 2000, 84, 5544 CrossRef CAS PubMed .
  79. D. J. Wales, Phys. Rev. E, 2017, 95, 030105(R) CrossRef PubMed .
  80. J. P. K. Doye and D. J. Wales, Phys. Rev. Lett., 1998, 80, 1357 CrossRef CAS .

Footnote

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

This journal is © The Royal Society of Chemistry 2018