The transformation mechanisms among cuboctahedra, Ino's decahedra and icosahedra structures of magic-size gold nanoclusters

Abstract

Gold nanoclusters possess multiple competing structural motifs with small energy differences, enabling structural coexistence and interconversion. Using a high-accuracy machine learned potential trained on some 20 000 density functional theory reference data points, we investigate transformation pathways connecting both high-symmetry and distorted cuboctahedra, Ino's decahedra and icosahedra for Au₅₅, Au₁₄₇, Au₃₀₉ and Au₅₆₁ nanoclusters. Our saddle point searches reveal that high-symmetry transformations from cuboctahedra and Ino's decahedra to icosahedra proceed through a single barrier and represent soft-mode-driven jitterbug-type and slip-dislocation motions. In addition, we identify lower-barrier asymmetric transformation pathways that drive the system into disordered, Jahn–Teller–stabilized distorted icosahedra. Minima Hopping sampling further uncovers, in this context, many such low-symmetry minima. Some of the newly identified global minima for Au₃₀₉ and Au₅₆₁ have energies that are up to 2.8 eV lower than the previously reported global minima. Hence, both the shapes and the transformation pathways studied in previous investigations are not the physically relevant ones. In contrast to the previously studied pathways, our transformation pathways give reasonable transformation times that are in rough agreement with experiments.

---

1 Introduction

Gold nanoclusters are of great interest due to their wide range of properties, enabling applications in many fields such as catalysis and biosensing.^1–3 In addition to the geometric ground state, gold nanoclusters can also adopt a large number of low-energy metastable structures with properties that are possibly quite different from those of the ground state. These meta-stable structures can either be defective structures of the underlying ground-state structural motif or be based on completely different structures.⁴ Theoretically, for instance, it was found that a truncated octahedron Au₂₀₁ is only 0.007 eV higher in energy than an icosahedron-like structure with five-fold symmetry.⁵ These small differences in energy between different structures are broadly consistent with experimental observations by the Palmer group.⁶

Larger nanoclusters are particularly stable if they consist of filled geometric shells. Such sizes are then called magic sizes. Icosahedra (I_h) with 55, 147, 309, and 561 atoms consist of two, three, four, and five shells around the central atom and are thus magic sizes. For the same number of atoms, one can also obtain Ino's decahedra (I-D_5h) and cuboctahedra (O_h) with filled shells.

Foster et al.⁶ and Koga et al.⁷ concluded from their experiments that for Au₅₆₁ nanoclusters, both the I-D_5h and face-centered cubic (FCC) structures, represented by the O_h motif, are energetically more favorable than I_h, with I-D_5h being only 0.040 eV higher in energy than the FCC. Furthermore, they observed that at temperatures above 300 K, O_h and I-D_5h motifs are more prevalent than I_h on carbon supports. They also observed amazing dynamic inter-conversions between I_h, I-D_5h and FCC-like structures occur on the time scale of seconds for nanoclusters with different numbers of atoms, such as Au₃₀₉ (ref. 8) and Au₅₆₁,⁶ highlighting the fluxional nature of these systems under experimental conditions.

In contrast, Density Functional Theory (DFT) calculations⁹ on isolated, ordered gold nanoclusters show that among the I_h, I-D_5h, and O_h structures, the I_h motif has the lowest energy, followed by I-D_5h and O_h. This trend is attributed to the presence of high-energy (100) facets in O_h and I-D_5h, making them energetically unfavorable, while rendering I_h the more stable structure.^10,11 It is noteworthy that, although the I_h motif represents the ground state for magic-size Lennard-Jones (LJ) clusters and is the lowest-energy configuration among ordered structures, it does not correspond to the true ground state of gold nanoclusters. Instead, the Global Minima (GM) are structures closely related to the I_h geometry but incorporating slight distortions. This behavior was first identified for the Au₅₅ nanocluster¹² and subsequently observed for larger nanoclusters, including Au₁₄₇, Au₃₀₉, and Au₅₆₁.⁴

The nature of the Au–C interaction explains the experimentally observed structural preferences. Gold on carbon supports typically exhibits weak binding energies of the order of ≈30–50 meV per atom, characteristic of physisorption with equilibrium separations around 3 Å.^13,14 Under these conditions, both geometric constraints imposed by the carbon lattice and interfacial interactions play a decisive role in determining nanocluster morphology. Since the cluster structures are almost degenerate in energy, such a relatively weak binding to a substrate can still affect the energy ordering of the nanoclusters.

In particular, graphene and graphite substrates favour the formation of close-packed FCC or mixed FCC–HCP structures, while suppressing non-crystalline motifs such as I_h. This behavior originates from the geometric matching between the hexagonal carbon lattice and the triangular arrangement of close-packed Au(111) layers, which enables flatter interfaces and better registry with the substrate. Consequently, supported nanoclusters tend to adopt layered geometries with predominantly (111) facets in contact with the carbon surface, while higher-energy facets such as (100) become unfavorable at the interface.¹⁵

Given the fact that I_h, I-D_5h, and O_h are fundamentally different, one might expect that the transformation from one structural motif to another proceeds by a nucleation process in which a large number of intermediate states are visited during the transformation. Schebarchov et al. calculated transformation pathways within Au₅₅, Au₈₅, and Au₁₄₇ nanoclusters.¹¹ Their paths cross several barriers while visiting about a dozen intermediate states with the lowest overall barriers of about 0.5 eV for Au₅₅. They interpret this transformation pathway as partial melting followed by crystallization.

Plessow investigated the atomistic transformation pathways from O_h to I_h structures in copper and other metal nanoclusters and discovered that the transition can proceed through the so-called jitterbug transformation, requiring only a single barrier.¹⁶ He showed that for gold nanoclusters smaller than Au₃₀₉, this barrier essentially vanishes, while for the Au₅₆₁, it is approximately 2.8 eV at the DFT PBE + D3 level, based on a saddle point obtained with the Gupta potential.¹⁷

In this work, the meta-stable and ground state structures of gold nanoclusters with the magic sizes of 55, 147, 309, and 561 are re-examined. Then the transformations of O_h, I-D_5h, and I_h gold nanoclusters into each other, referred to as high-symmetry transformations from now on, are studied. We consider not only transformations leading to the I_h structure but also those resulting in various distorted forms—namely, distorted I_h (d-I_h), distorted I-D_5h (d-I-D_5h), and distorted O_h (d-O_h). These will be referred to as asymmetric transformations throughout the remainder of this paper. Notably, similar distorted structures have previously been described as “amorphous” in ref. 4, although they retain identifiable symmetry motifs and are therefore characterized here as distorted.

2 Methods

DFT remains the method of choice for computing accurate energies and forces. Single-point calculations and geometry relaxations are now possible for nanoclusters containing hundreds of atoms. However, simulations requiring numerous energy and force evaluations—such as vibrational mode calculations, saddle point searches, or long-time Molecular Dynamics (MD) simulations remain computationally prohibitive for large nanoclusters on the DFT level. For this reason, such theoretical investigations of large nanoclusters were done nearly exclusively by approximate interatomic potentials such as the Gupta potential or MEAM^18,19 methods. The price for the gain in speed is reduced accuracy. Fig. 1 illustrates the poor agreement between the energies obtained from DFT on the one hand and either the Gupta or MEAM potential on the other hand over a diverse dataset of 19 383 gold nanoclusters.

Fig. 1 Correlation of the energies calculated by three different force fields for the structures in the training dataset of the MLIP.

The advances in the development of machine-learned interatomic potentials (MLIPs) have changed the situation dramatically. MLIPs can nowadays give DFT accuracy at a small fraction of the DFT cost, enabling complex simulations of large systems. We therefore use such an MLIP, namely NequIP,^20–22 to investigate large gold nanoclusters. The MLIP was trained on a dataset of 19 383 nanoclusters, each containing between 10 and 90 atoms. The DFT reference energies were computed with the Projector-Augmented Wave (PAW) method^23,24 using 11 valence electrons as implemented in the Vienna Ab initio Simulation Package (VASP)^25,26 together with the Perdew–Burke–Ernzerhof (PBE) exchange-correlation functional.²⁷ To remove interactions between periodic images, a vacuum of 12 Å was introduced in all directions around the nanoclusters.

As shown in Fig. 1, the resulting MLIP exhibited a root mean square error of 7.9 meV per atom for energies and 75.4 meV Å⁻¹ for forces, demonstrating excellent accuracy for the subsequent tasks. To further assess the transferability of the MLIP, we evaluated its performance on larger nanoclusters—Au₅₅, Au₁₄₇, Au₃₀₉, and Au₅₆₁—that were not included in the training set. The model yielded a root mean square error of 5.2 meV per atom for energies (see Fig. S3) and 37.9 meV Å⁻¹ for forces. MLIP forces and energies were used for MD, saddle point searches with both the Nudged Elastic Band (NEB)²⁸ and COMPASS,²⁹ Minima Hopping method (MH) runs^30–34 and vibrational mode analysis within the Atomic Simulation Environment (ASE).³⁵ The training set was generated iteratively using MD and MH simulations combined with active learning, with the MLIP training set expanded after each iteration. The resulting dataset spans an energy range of up to 1.12 eV per atom. The energy and size distributions of the training structures are shown in the SI, Fig. S1 and S2, respectively.

MH explores the PES by going over low barriers and visiting low-energy local minima in the pursuit of finding the GM. By steering the MH search towards a target structure, a collection of minima can be obtained to connect the initial structure to the target structure with a physically realizable pathway.³⁶ The steering is obtained by performing MH on a biased PES (biased MH) Ẽ(R) = E(R) + ωD(R, R^tar) where ω is the bias weight and D(R, R^tar) = |F(R) − F(R^tar)|² is a distance function that is zero if the current structure R is equal to the target structure R^tar and positive for all other structures.

The fingerprint vector F(R) contains the eigenvalues of the Laplace matrix L. Position-dependent matrix elements are given by:

Since the length of a fingerprint vector F is N for a nanocluster containing N atoms, it can actually not uniquely specify the 3N degrees of freedom of all the atoms. However, for a biasing scheme, this is not problematic.

To determine the transformation pathway between different atomic configurations of gold nanoclusters, NEB and COMPASS calculations were conducted. The mapping of atomic indices of the initial and final structures in COMPASS calculations for asymmetric transformations was determined by a simple Root Mean Square Distance alignment (RMSD).³⁷ The parameters of COMPASS are provided in table SV. For the high-symmetry transformations, a total of 21 images was picked, including the initial and final configurations for each NEB pathway from biased MH. The NEB calculations converged to a force tolerance of five meV Å⁻¹, using the Fast Inertial Relaxation Engine (FIRE),³⁸ ensuring high precision in the computed transition states and reaction pathways. To calculate the vibrational frequencies and modes, a five-point finite difference method with a finite displacement of 0.01 Å was employed to compute the Hessian matrix.

The low energy barriers determine the dynamics of a system at ambient temperatures. Empirical findings suggest that MD trajectories starting from a minimum along a soft vibrational mode are more likely to cross low-energy barriers compared to those starting in other directions.³⁹ The softest mode on the PES corresponds to the lowest non-zero vibrational frequency and reveals the most transformation-prone direction in configuration space. By aligning the initial velocities along this mode and scaling their magnitude to correspond to a specific temperature, one can control both the energy input and its directional bias, effectively triggering rapid structural transformations. In a related context, Wang et al.⁴⁰ studied titanium oxide nanoparticles under negative pressure and found also that soft vibrational modes can trigger coherent structural phase transformations.

According to classical transition state theory, the reaction rate, k, for crossing a barrier is given by:

is the attempt frequency that can be considered as the rate with which the system tries to overcome a barrier at temperature T with numerical value of 6.25 × 10¹² (Hz) at T = 300 K. F^‡_S(β) is the free energy of the saddle point and F_A(β) is the free energy of the local minimum. The ‡ symbol indicates the usage of positive modes only when calculating the free energy. The six zero modes resulting from rotational and translational symmetries do not contribute. In our calculations, we approximate the free energy by the energy.

3 Results and discussion

3.1 Equilibrium structures

To thoroughly sample the PES of Au₅₅, Au₁₄₇, Au₃₀₉ and Au₅₆₁, we used MH combined with MLIP without any biasing. From the ensemble of the visited minima, we selected the twelve lowest-energy candidates. The selected candidates were further relaxed up to one meV Å⁻¹ using VASP with ISPIN = 2 to enable spin-polarized calculation, ISMEAR = 0 for Gaussian smearing with width SIGMA = 0.0005 eV. The small SIGMA value is chosen because we are interested in the zero temperature limit where Jahn–Teller distortions occur in our nanoclusters.⁴¹

Table 1 shows that the newly identified global minima of Au₃₀₉ (Fig. 2b) and Au₅₆₁ (Fig. 2e) are 0.260 and 2.896 eV lower in energy than the previously reported structures,^4,42 even though their overall morphologies remain largely unchanged. The structures have I_h character, displaying well-defined triangular facets on one side and a somewhat distorted character on the other side. We found some 200 structures of this d-I_h type, which suggests that the I_h funnel dominates the low-energy landscape. Alternative motifs remain, however, competitive in energy (Table 1). The coordinates of the 10 lowest-energy distorted structures are provided in the SI.

Table 1 Energy differences relative to the geometry optimized I_h structure (in eV). Our lowest energy structures are contrasted with those from ref. 4

Rank	Au55	Au147	Au309	Au561
a Ref. 4 does not present the structures of Au55 and Au147.
d-Ih 1	−2.387	−4.382	−8.625	−13.455
d-Ih 2	−2.385	−4.304	−8.574	−13.395
d-Ih 3	−2.353	−4.243	−8.412	−13.334
d-Ih 4	−2.318	−4.189	−8.396	−13.306
d-Ih 5	−2.311	−4.187	−8.224	−13.305
Ref. 4	^a	^a	−8.366	−10.559
Oh	1.085	3.110	4.802	3.536
I-D5h	0.434	2.518	1.173	1.024

---

Fig. 2 New GM structures for Au₃₀₉ (top row) and Au₅₆₁ (bottom row). (a and d) Atoms whose positions coincide with those of an I_h are shown in gold. (b and e) Complete nanocluster structures, where ordered (I_h-matching) atoms are shown in golden color and distorted atoms are shown in red. (c and f) Isolated views of the distorted (red) atoms extracted from the corresponding structures.

Fig. 3 shows that while a well-trained MLIP can capture the overall features of the energetic landscape of gold nanoclusters, it lacks quantum-mechanical principles to account for subtle electronic effects such as a lowering of the occupied levels or Jahn–Teller distortions, which are known to play a decisive role in stabilizing lower-energy isomers via symmetry breaking.⁴³ Such effects are, for instance, observed when the MLIP-relaxed I_h was tightly DFT post-relaxed while its electronic structure was calculated using VASP with ISPIN = 2, ISMEAR = 0 and SIGMA = 0.0005.

Fig. 3 Energy levels near the Fermi level calculated for three different structures of Au₅₅-I_h obtained with three approaches: (i) scaling LJ₅₅-I_h cluster on the MLIP PES with a factor (scaled LJ), (ii) relaxing the scaled LJ structure on the MLIP PES (MLIP-relaxed), and (iii) further relaxing the MLIP-relaxed structure on the DFT PES (VASP-relaxed). The scaling factor in approach (i) was the one that gave the lowest MLIP energy. Occupied and unoccupied levels are shown as solid and dashed lines, respectively, while the Fermi level is indicated by a black dotted line.

While the MLIP-relaxed I_h structure preserves degenerate electronic levels, these are lifted upon further relaxation of forces up to one meV Å⁻¹ in spin-polarized DFT, indicating a Jahn–Teller distortion.⁴³ This lifting of degeneracy enables complete filling of the resulting split levels, as shown in Fig. 3 for Au₅₅-I_h. The accompanying changes in the Au–Au bond lengths and the resulting broadening of their distribution are shown in Fig. 4. The total energy of the MLIP-relaxed Au₅₅-I_h nanocluster decreased by 0.046 eV during this Jahn–Teller distortion, and the largest atomic displacement observed during this step was 0.053 Å. The breaking of the full I_h symmetry is therefore hard to detect by eye, and we will continue to refer to the Jahn–Teller distorted structure as the I_h for the rest of the paper.

Fig. 4 Histogram comparing the bond length distributions of the Au₅₅-I_h nanocluster as obtained from scaled LJ, MLIP-relaxed, and VASP-relaxed structures.

The same effect exists for all magic I_h sizes, although the degree of gap opening decreases with size as the density of states near the Fermi level increases. So, the rule established for smaller nanoclusters,⁴⁴ namely that nanoclusters adopt a structure that allows them to fill degenerate levels completely, is also valid for larger nanoclusters, even though their overall shape is dictated by geometric effects.

In analogy to the I_h, where distorted I_h nanoclusters were lower in energy than the I_h, there also exist d-I-D_5h and d-O_h that are lower in energy than their high symmetry structures. Four such structures are shown in Fig. 5. The number of distorted I-D_5h and O_h that are lower in energy than their high-symmetry counterparts is smaller than in the case of the I_h.

Fig. 5 Examples of two distorted O_h-Au₁₄₇ (a and b) and two I-D_5h-Au₁₄₇ (c and d). Red spheres indicate atoms that deviate from their ideal positions in the corresponding high-symmetry nanocluster.

3.2 Minimum energy transformation pathways

Performing biased MH runs where the target structure was I_h, revealed immediately the so-called Jitterbug transformation for elastic bodies connecting O_h → I_h and a cooperative slip-dislocation mechanism reported for an elastic body by Koga et al. for the I-D_5h → I_h transformation.⁷ In this way, we could show that these transformations, which were only defined for shapes, can also be represented by simple and short atomic trajectories.

To assess the accuracy of the MLIP for structural transformations, NEB calculations were performed for Au₅₅ and Au₁₄₇ on the DFT PES with an initial guess consisting of nine images selected from biased MH. For comparison, NEB calculations were also performed with 21 images on the MLIP PES. The NEB was relaxed up to ten meV per atom on DFT PES and five meV per atom on MLIP PES. Solid lines in Fig. 6 show the energy of NEB images on the DFT PES with no climbing image, i.e., intermediate image that feels no forces from the NEB springs. Dash-dotted lines represent the MLIP energy profiles for Au₅₅ and Au₁₄₇ nanoclusters along high-symmetry transformations, including a climbing image. The agreement between the two methods confirms that the MLIP reproduces DFT-level energetics throughout these transformations within the MLIP error range of 7.9 meV per atom. Having established the accuracy of MLIP, we subsequently employed the MLIP to perform NEB calculations using 21 equally spaced structures for high-symmetry transformations of Au₃₀₉ and Au₅₆₁ nanoclusters.

Fig. 6 Comparison of energy profiles along the transformation pathways for Au₅₅ and Au₁₄₇ nanoclusters obtained from NEB calculations for O_h → I_h and I-D_5h → I_h transformations (top and bottom, respectively). The energy of structures along the transformation pathway remains nearly the same, independently of whether single-point DFT energies of MLIP-relaxed images are calculated (dashed line) or the NEB is fully relaxed directly at the DFT level (solid line).

Fig. 7 shows that NEB paths indeed connect the two structures while crossing only one barrier. However, initial vibrational calculations using MLIP indicated that the identified saddle points, i.e., the climbing images of the NEB, were not first-order. In this case, it should be possible to find lower-energy first-order saddle points. Such lower saddle points on the high-symmetry transformation pathway were, however, not found. To eliminate the possibility that the second-order character of the saddle points is an artifact of numerical noise, which is always present in MLIP or DFT calculations, and to isolate purely geometric effects, we additionally performed NEB simulations using a simple LJ potential, which is noise-free.

Fig. 7 The NEB pathway calculated by MLIP, dash-dotted line, for O_h → I_h and I-D_5h → I_h transformation, top and bottom, respectively. The single-point energies of NEB individual images calculated on DFT PES are represented by a dotted line. In the insets, the relevant NEB energy pathway for geometrically scaled and relaxed structures on the LJ potential is presented.

The tests were carried out on scaled geometries of the gold I-D_5h, O_h, and I_h structures. The LJ transformation pathways were virtually (up to a scaling factor) identical to those of the gold nanoclusters, showing that the jitterbug and slip-dislocation transformations are purely geometrical and largely independent of the interaction potential. A similar behavior is observed for ionized Au₅₅: the geometric pathway of the high-symmetry transformation remains essentially unchanged, while the barrier heights remain small (see Section S3). In the limit of large nanoclusters, the transformation pathways and barrier heights are expected to become even less sensitive to the charging of the nanocluster. The additional low-curvature directions with multiple negative eigenvalues, observed on the MLIP PES, were not present when the Hessian was calculated analytically for the corresponding LJ saddle points. This strongly suggests that the extra negative modes in the MLIP results stem from numerical noise in the MLIP PES.

It is also worth noting in this context that displacements along the lowest positive curvature direction of the LJ saddle points revealed a wide and nearly flat valley in the saddle point region with a shape similar to the one of the gold PES depicted in Fig. 8 and 9. This flatness can introduce significant numerical noise when estimating vibrational frequencies via finite differences on the MLIP PES, complicating the accurate characterization of saddle points. The energy profile along this direction, computed for the saddle points of the LJ and gold nanoclusters, is provided in Fig. S6.

Fig. 8 Energy for displacements along vibrational modes of Au₅₅-I-D_5h. The solid lines show the energy profiles when displacing atoms along the ten non-zero lowest-frequency (softest) modes. For comparison, the dotted lines of the same color indicate the energy of the corresponding harmonic potential for each mode. The dash-dotted blue line represents the energetic distance of the highest-frequency vibrational modes, ħω_max. The fact that this line is at room temperature below k_BT ≈ 0.025 eV shows that all vibrational modes are thermally accessible at room temperature.

Fig. 9 Energy for displacements along vibrational modes of the saddle point of the Au₅₅-I-D_5h → Au₅₅-I_h transformation. The conventions are the same as in Fig. 8. The additional dashed line denotes the imaginary mode at the saddle point, arising from a negative curvature of the PES along this mode.

As can be seen from Fig. 7, the barrier height of small nanoclusters is almost negligible, but it increases as the nanoclusters grow larger, reaching approximately 4.2 eV for Au₅₆₁-O_h. The relaxed pathways indicate that the transformations from O_h and I-D_5h to I_h are strongly exothermic, with the degree of exothermicity decreasing as the nanocluster size increases. In our calculations, the energy pathway for the O_h → I_h transformation (Fig. 7) remains qualitatively consistent across all investigated nanocluster sizes. The energy profiles obtained here are similar to the transformation pathway reported by Plessow¹⁶ for copper nanoclusters, although the energy barriers we observe are higher. Specifically, our calculated activation energies for the O_h → I_h transformation are higher than those reported by Plessow by 0.440, 1.698, and 1.402 eV for Au₁₄₇, Au₃₀₉, and Au₅₆₁, respectively. This difference can most likely be attributed to the improved accuracy of the MLIP employed in this work compared to the Gupta potential used by Plessow.¹⁶

The NEB pathway of the LJ system, shown in the insets of Fig. 7, also follows the same trend as the gold MLIP results despite the absence of any bond-bending terms in the pairwise LJ potential. This indicates that the transformations are predominantly geometry-driven and largely independent of the material.

Fig. 10 provides a visual representation of the minimum energy pathway for the high-symmetry transformations of Au₅₅ alongside the lowest vibrational modes. The results for other sizes, shown in Fig. S8, follow almost the same trend as Au₅₅ nanoclusters. The xyz-coordinates of atoms for the transformation pathways are provided in the SI. In Fig. 10, one can clearly see that the high-symmetry transformation is not a nucleation process. Instead, it is a concerted, smooth, simultaneous movement of all atoms. At the starting local minimum, the tangent to the minimum energy pathway is well aligned with the lowest-frequency vibrational mode, indicating that the system predominantly evolves along the softest mode direction during a high-symmetry transformation.

Fig. 10 (a) Depicts the high-symmetry Au₅₅-I-D_5h → Au₅₅-I_h transformation pathway (side view), while (b) illustrates the high-symmetry Au₅₅-O_h → Au₅₅-I_h pathway (viewed along (100)). Golden spheres represent the initial structures (I-D_5h or O_h), green spheres indicate the final structure (I_h), red lines show the direction of the softest mode, blue lines trace the NEB transformation pathway, and orange highlights the saddle-point configuration. The xyz-coordinates of these pathways are provided in the SI.

To determine the lowest barrier that a high-symmetry nanocluster must overcome to transform into the GM, we employed the COMPASS method.²⁹ The complex transformation pathway identified by COMPASS clarified why the NEB method fails when applied to the asymmetric transformation. Interestingly, barriers leading into these low-energy distorted structures were lower than those leading to the I_h. As a result, it is much more likely that the system will transform into a d-I_h than into the I_h. Experimentally, one can presumably not distinguish the two forms of I_h.

Pathways that lead into the GM are complicated, as shown in Fig. 11. Multiple barriers of various heights have to be crossed to reach the GM. To quantify the symmetry breaking associated with the stabilization of these d-I_h structures relative to the I_h geometry, we analyze the energy and the structural deviation from the I_h, measured by RMSD from the I_h structure. As shown in Fig. 12, the energy decreases with increasing RMSD, indicating that lower-energy configurations are obtained via distortions of the I_h structure. The magnitude of these distortions remains limited, such that the resulting configurations preserve the overall I_h-like morphology while lacking exact symmetry.

Fig. 11 Energy profiles along the COMPASS transformation pathways are shown as colored solid lines. The top, middle, and bottom panels correspond to pathways starting from O_h, I-D_5h, and I_h, respectively. Each curve displays the energy difference relative to the starting structure of its panel. The dotted black line denotes zero energy. The dash-dotted magenta line marks the energy of the highest-energy structure along the pathway, indicating the upper bound of the barrier. The dashed black line shows the catchment basin visited during the transformation, obtained by relaxing the COMPASS pathway structures on the MLIP PES. Colors on the solid line indicate the structural motif: red for O_h, green for I-D_5h, blue for I_h, and orange for fully distorted nanoclusters.

Fig. 12 Comparison of energy difference relative to the I_h and RMSD from the I_h structure for Au₅₅ nanoclusters. Points represent minima sampled within the I_h funnel. The GM (d-I_h) corresponds to a d-I_h structure with lower energy than the I_h.

The asymmetric transformations do not proceed along a single dominant vibrational mode. When expressing the initial atomic displacement by a linear combination of the vibrational modes, soft modes contribute most significantly to this motion, whereas higher-frequency (hard) modes play a minor role. To quantify this relationship, we compute the angle between the vector connecting the initial structure to the asymmetric transformation saddle point and each vibrational eigenmode. The resulting angular distribution, shown in Fig. S7 illustrates how individual vibrational modes participate in shaping the reaction direction at the beginning of asymmetric transformations.

In Table 2, the movement amplitudes of atoms during various transformations are presented. In the inner shells, the atoms move only by a small fraction of a bond length, while in the outer shells, the movement is larger, reaching 1.15 times the gold bond length in nanoclusters, i.e., 2.94 Å.

Table 2 Minimum and maximum atomic displacements (in Å) observed during various transformations, excluding the innermost atom

Transformation		Au55	Au147	Au309	Au561
Oh → Ih	Min	0.4	0.4	0.4	0.4
	Max	1.3	1.9	2.6	3.2
Oh → GM	Min	0.1	0.1	0.1	0.1
	Max	2.6	2.4	2.7	3.2
I-D5h → Ih	Min	0.1	0.2	0.2	0.2
	Max	1.5	2.3	3.1	3.9
I-D5h → GM	Min	0.1	0.1	0.1	0.1
	Max	2.3	2.9	2.9	3.0
Ih → GM	Min	0.1	0.1	0.1	0.1
	Max	2.8	3.1	3.3	3.4

---

Portales et al. attributed a ∼200 GHz shift in low-frequency Raman scattering to quadrupolar (E_g and T_2g) vibrational modes in 4.3–5.3 nm gold nanocrystals.⁴⁵ Although these nanocrystals contain far more atoms than our nanoclusters, their reported T_2g mode aligns with the softest vibrational modes of the O_h structures listed in Table 3. Notably, the T_2g mode described by Portales et al. corresponds closely to the jitterbug motion associated with the O_h → I_h transformation. Still, it differs significantly from the slip-dislocation twisting involved in the I-D_5h → I_h transformation. The red lines, direction of O_h softest mode, in Fig. 10b align with the arrows showing the T_2g mode in the Portales paper.

Table 3 Energy (in meV) and frequencies (in GHz) of the softest modes for O_h, I-D_5h and I_h

	Oh		I-D5h		Ih
	E	f	E	f	E	f
Au55	0.8	200.86	1.2	290.80	2.5	602.58
Au147	1.2	287.80	1.1	263.81	2.1	500.65
Au309	1.1	278.80	0.9	227.84	1.7	398.72
Au561	1.1	254.82	0.8	200.86	1.3	323.77

---

3.3 Kinetics of transformation

The barriers for all kinds of transformations in the small Au₅₅ and Au₁₄₇ nanoclusters are quite low, as shown in Table 4. Classical transition state theory gives transformation times of the order of picoseconds (ps) for these barriers, consistent with predictions by Schebarchov et al.¹¹ Experimental transformation rates are, however, of the order of seconds. This discrepancy is likely due to the interaction between the nanoclusters and the substrate, as Foster et al. noted in their SI.⁶ Presumably, the interaction of the nanocluster with the substrate makes rearrangements harder. For the larger nanoclusters, the transition times explode, whereas the experimental time scales remain nearly constant.⁴⁶ It is hard to believe that the interaction with a substrate can speed up transformations in a nanocluster that is unwilling to transform.

Table 4 Barrier heights (in eV) for both high-symmetry and asymmetric transformations. The shows the barrier heights extracted from Van't Hoff plots based on the standard MD data (see Fig. 13). Values in square brackets represent the resulting lifetime of initial nanoclusters in ps using transition-state theory at 300 K for these barriers

Transformation	Au55	Au147	Au309	Au561
a Au55-Oh transforms into Au55-Ih instead of d-Oh.
Oh → Ih	2.09 × 10^−2	4.40 × 10^−1	1.70 × 10^0	4.20 × 10^0
	3.6 × 10^−1	3.9 × 10^6	5.8 × 10^27	5.8 × 10^69
Oh → d-Oh	2.09 × 10^−2^a	2.28 × 10^−1	2.91 × 10^−1	3.62 × 10^−1
	3.6 × 10^−1	1.1 × 10^3	1.2 × 10^4	1.9 × 10^5
	4.55 × 10^−2	3.33 × 10^−1	2.44 × 10^−1	1.74 × 10^−1
	9.3 × 10^−1	6.3 × 10^4	2.0 × 10^3	1.3 × 10^2
I-D5h → Ih	1.19 × 10^−1	4.84 × 10^−1	1.36 × 10^0	2.95 × 10^0
	1.6 × 10^1	2.1 × 10^7	1.1 × 10^22	5.8 × 10^48
I-D5h → d-D5h	9.96 × 10^−2	1.60 × 10^−1	1.92 × 10^−1	2.15 × 10^−1
	7.6 × 10^0	7.8 × 10^1	2.7 × 10^2	6.6 × 10^2
	1.74 × 10^−1	2.29 × 10^−1	2.29 × 10^−1	2.13 × 10^−1
	1.4 × 10^2	1.1 × 10^3	1.1 × 10^3	6.2 × 10^2
Ih → Oh	1.55 × 10^0	4.23 × 10^0	8.26 × 10^0	13.22 × 10^0
	2.2 × 10^25	1.8 × 10^70	9.2 × 10^137	2.0 × 10^221
Ih → I-D5h	1.08 × 10^0	2.99 × 10^0	5.68 × 10^0	8.82 × 10^0
	2.4 × 10^17	3.8 × 10^49	5.3 × 10^94	3.3 × 10^147
Ih → d-Ih	4.13 × 10^−2	1.75 × 10^−1	2.06 × 10^−1	5.44 × 10^−1
	7.9 × 10^−1	1.4 × 10^2	4.7 × 10^2	2.2 × 10^8
	4.36 × 10^−2	3.08 × 10^−1	3.68 × 10^−1	3.59 × 10^−1
	8.6 × 10^−1	2.4 × 10^4	2.4 × 10^5	1.7 × 10^5

---

How can this be explained? One possibility would be that the free energy differences between the saddle point and the initial minimum are at room temperature much smaller than the energy difference at zero temperature that we have used in our calculations. Since soft modes lower the free energy more than hard modes, we have studied the modes at the Au₅₅-I-D_5h minimum and its saddle point toward the Au₅₅-I_h. In Fig. 8 and 9, we plot both the exact potential and its quadratic approximation along the vibrational modes. Atomic displacements are generated according to

R^displaced_i = R⁰_i + xν_i, (3)

where R⁰_i denotes the initial atomic positions of atom i, ν is an vibrational mode vector and x is the displacement amplitude since

For each mode, atoms are displaced parallel and anti-parallel to the mode by an amplitude x to map out the local energy landscape around the reference structure. To assess whether the displacement amplitudes remain within the validity region of the harmonic approximation, the harmonic energy profile along the corresponding mode, given by , is shown as a reference using dotted curves. m is the mass of a gold atom. By comparing Fig. 8 and 9, one sees that the saddle point has somewhat softer modes than the minimum, but the effect is too weak to bring down the theoretical predictions to the experimental values.

The second possibility is that the systems in the experiment do not follow the symmetric jitterbug and slip-dislocation transformations but undergo asymmetric transformations, which have, as we have seen, a much lower barrier. To check this hypothesis, we performed ordinary MD trajectories starting from O_h, I-D_5h, and I_h. The trajectories were interrupted by local geometry relaxations every 50 time steps, i.e., every 250 fs, to check if the trajectory had left its initial catchment basin. If this were the case, the MD was stopped, and the time was recorded. Performing this procedure at different temperatures provided us with the data necessary to make the Van't Hoff plot shown in Fig. 13.

Fig. 13 Van't Hoff plots for structural transformations of Au₅₅, Au₁₄₇, Au₃₀₉ and Au₅₆₁ starting from O_h (top), I-D_5h (middle) and I_h (bottom). The nearly identical slopes across different nanocluster sizes indicate that the transformation barriers are effectively size-independent.

For small nanoclusters such as Au₅₅ and Au₁₄₇, the barrier heights obtained from the Van't Hoff analysis show good agreement with the calculated barriers along high-symmetry transformation pathways. This consistency indicates that the transformation kinetics in this size range are close to the high-symmetry transformations and well described by classical transition state theory.

For large nanoclusters, eqn (2) together with the barrier heights from Table 4, predicts transformation times for high-symmetry transformations that are longer than the age of the universe. Fig. 13 shows, however, that the effective barrier heights extracted from the Van't Hoff plots are nearly independent of nanocluster size, which is in agreement with experiments and COMPASS calculations. This implies that the larger nanoclusters never follow the symmetric jitterbug and slip-dislocation transformations. Instead, they follow asymmetric transformation pathways that consist of a series of localized transformation steps over low-energy barriers. One such step is shown in Fig. 15 for all nanocluster sizes. The localized nature becomes clearly visible by comparing the atomic displacements of these transformations with the ones shown in Fig. 15 and S8.

The size-independent nature of these asymmetric rearrangements is further illustrated in Fig. 14. Across Au₅₅, Au₁₄₇, Au₃₀₉, and Au₅₆₁, the distribution of barrier heights remains remarkably similar. The large number of nearly degenerate d-I_h structures separated by low barriers indicates that large regions in the configuration space are accessible in the I_h funnel. Since the transformation times are, in most cases, much shorter than the time resolution of the experimental apparatus⁶ used to investigate the structure of these nanoclusters, experimentally measured properties will, in general, be averages over several distorted structures.

Fig. 14 Barrier heights for the d-I_h → d-I_h transformation for Au₅₅, Au₁₄₇, Au₃₀₉, and Au₅₆₁.

Fig. 15 Localized displacement of atoms when transferring from one d-I_h structure to another d-I_h for Au₅₅ (a), Au₁₄₇ (b), Au₃₀₉ (c) and Au₅₆₁ (d).

Conclusions

In this work, we first fully mapped out the transformation pathways between I_h, I-D_5h, and O_h gold nanoclusters with 55, 147, 309, and 561 atoms using a highly accurate machine learned interatomic potential. We find that these transformations are concerted movements of all atoms and that they proceed along the softest vibrational modes of the nanocluster. The atoms travel only relatively small distances during these transformations. The transformations are the atomistic analogues of the jitterbug and slip-dislocation mechanisms of elastic bodies. The height of the barrier of these transformations increases with increasing system size to values that make the transformation impossible on any experimental time scale. Hence, these symmetric transformation pathways cannot be the relevant ones in nature. We found numerous novel d-I_h structures that are much lower in energy than the I_h, as well as low-lying distorted I-D_5h and O_h structures. We found that the barriers of transformations among these distorted structures are much lower than the barriers between the high-symmetry structures. This is because these transformations involve typically only a subset of atoms within a small region of the nanocluster. As a consequence, they can be overcome on relatively short time scales. We therefore postulate that the experimentally observed transitions between I_h, I-D_5h, and O_h shapes are not transitions between the high-symmetry shapes but among distorted-type structures.

Conflicts of interest

There are no conflicts to declare.

Data availability

The coordinate files for all structures and transformations reported in this article have been deposited and are provided in the accompanying zipped archive. Further information regarding the data can be obtained by contacting the corresponding author via email.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d6na00012f.

Acknowledgements

We thank Prof. Richard Palmer for the interesting discussions. Financial support for this project was provided by SNF under grant number 200021_191994. This work was supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project ID lp08 on Alps. Calculations were also performed at sciCORE (https://scicore.unibas.ch/), the scientific computing center at the University of Basel.

Notes and references

1. M.-C. Daniel and D. Astruc, Chem. Rev., 2004, 104, 293–346.
2. P. L. Rodríguez-Kessler and A. Muñoz-Castro, Inorg. Chim. Acta, 2025, 574, 122376.
3. M. Vidotti, R. F. Carvalhal, R. K. Mendes, D. Ferreira and L. T. Kubota, J. Braz. Chem. Soc., 2011, 22, 3–20.
4. S. Jindal and S. S. Bulusu, J. Chem. Phys., 2020, 152, 154302.
5. K. Bao, S. Goedecker, K. Koga, F. Lançon and A. Neelov, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 041405.
6. D. Foster, R. Ferrando and R. Palmer, Nat. Commun., 2018, 9, 1323.
7. K. Koga, T. Ikeshoji and K.-i. Sugawara, Phys. Rev. Lett., 2004, 92, 115507.
8. M. Dearg, C. Roncaglia, D. Nelli, E. Y. El Koraychy, R. Ferrando, T. J. A. Slater and R. E. Palmer, Nanoscale Horiz., 2024, 9, 143–147.
9. H. Li, L. Li, A. Pedersen, Y. Gao, N. Khetrapal, H. Jónsson and X. C. Zeng, Nano Lett., 2015, 15, 682–688.
10. F. Baletto and R. Ferrando, Rev. Mod. Phys., 2005, 77, 371–423.
11. D. Schebarchov, F. Baletto and D. Wales, Nanoscale, 2018, 10, 2004–2016.
12. 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–1603.
13. G. Giovannetti, P. A. Khomyakov, G. Brocks, V. v. Karpan, J. van den Brink and P. J. Kelly, Phys. Rev. Lett., 2008, 101, 026803.
14. C. Gong, G. Lee, B. Shan, E. M. Vogel, R. M. Wallace and K. Cho, J. Appl. Phys., 2010, 108, 123711.
15. M. Settem, M. M. Gianetti, R. Guerra, N. Manini, R. Ferrando and A. Giacomello, Small Sci., 2024, 4, 2400078.
16. P. N. Plessow, Phys. Chem. Chem. Phys., 2020, 22, 12939–12945.
17. R. P. Gupta, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 23, 6265.
18. M. S. Daw and M. I. Baskes, PPhys. Rev. B: Condens. Matter Mater. Phys., 1984, 29, 6443.
19. C. O'Brien, C. Barr, P. Price, K. Hattar and S. Foiles, J. Mater. Sci., 2018, 53, 2911–2927.
20. S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt and B. Kozinsky, Nat. Commun., 2022, 13, 2453.
21. B. Kozinsky, A. Musaelian, A. Johansson and S. Batzner, Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, 2023, pp. 1–12.
22. C. W. Tan, M. L. Descoteaux, M. Kotak, G. d. M. Nascimento, S. R. Kavanagh, L. Zichi, M. Wang, A. Saluja, Y. R. Hu and T. Smidt, Digital Discovery, 2026, 1558–1567.
23. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953.
24. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
25. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
26. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
27. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
28. H. Jónsson, G. Mills and K. W. Jacobsen, Classical and Quantum Dynamics in Condensed Phase Simulations, World Scientific, 1998, pp. 385–404.
29. M. Sommer-Jörgensen and S. Goedecker, J. Chem. Phys., 2024, 160, 064110.
30. S. E. Schönborn, S. Goedecker, S. Roy and A. R. Oganov, J. Chem. Phys., 2009, 130, 144108.
31. S. Goedecker, J. Chem. Phys., 2004, 120, 9911–9917.
32. S. Roy, S. Goedecker and V. Hellmann, Phys. Rev. E, 2008, 77, 056707.
33. M. Amsler and S. Goedecker, J. Chem. Phys., 2010, 133, 224104.
34. M. Krummenacher, M. Gubler, J. A. Finkler, H. Huber, M. Sommer-Jörgensen and S. Goedecker, SoftwareX, 2024, 25, 101632.
35. A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer and C. Hargus, et al., J. Phys.: Condens. Matter, 2017, 29, 273002.
36. D. S. De, M. Krummenacher, B. Schaefer and S. Goedecker, Phys. Rev. Lett., 2019, 123, 206102.
37. M. Griffiths, S. P. Niblett and D. J. Wales, J. Chem. Theor. Comput., 2017, 13, 4914–4931.
38. E. Bitzek, P. Koskinen, F. Gähler, M. Moseler and P. Gumbsch, Phys. Rev. Lett., 2006, 97, 170201.
39. M. Sicher, S. Mohr and S. Goedecker, J. Chem. Phys., 2011, 134, 044106.
40. K. Wang, C. Molteni and P. D. Haynes, Nano Lett., 2022, 22, 5922–5928.
41. I. B. Bersuker, The Jahn–Teller Effect, Cambridge Cambridge University Press, 2006.
42. I. Garzón, K. Michaelian, M. Beltrán, A. Posada-Amarillas, P. Ordejón, E. Artacho, D. Sánchez-Portal and J. Soler, Phys. Rev. Lett., 1998, 81, 1600.
43. H. A. Jahn and E. Teller, Proc. Roy. Soc. Lond. Math. Phys. Sci., 1937, 161, 220–235.
44. G. Fisicaro, B. Schaefer, J. A. Finkler and S. Goedecker, Mater. Adv., 2023, 4, 1746–1768.
45. H. Portales, N. Goubet, L. Saviot, S. Adichtchev, D. Murray, A. Mermet, E. Duval and M.-P. Pileni, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 14784–14789.
46. Z. Wang and R. Palmer, Phys. Rev. Lett., 2012, 108, 245502.

---