Impurity effects on solid–solid transitions in atomic clusters

Abstract

We use the harmonic superposition approach to examine how a single atom substitution affects low-temperature anomalies in the vibrational heat capacity (C_V) of model nanoclusters. Each anomaly is linked to competing solidlike “phases”, where crossover of the corresponding free energies defines a solid–solid transition temperature (T_s). For selected Lennard-Jones clusters we show that T_s and the corresponding C_V peak can be tuned over a wide range by varying the relative atomic size and binding strength of the impurity, but excessive atom-size mismatch can destroy a transition and may produce another. In some tunable cases we find up to two additional C_V peaks emerging below T_s, signalling one- or two-step delocalisation of the impurity within the ground-state geometry. Results for Ni₇₄X and Au₅₄X clusters (X = Au, Ag, Al, Cu, Ni, Pd, Pt, Pb), modelled by the many-body Gupta potential, further corroborate the possibility of tuning, engineering, and suppressing finite-system analogues of a solid–solid transition in nanoalloys.

---

I. Introduction

Multi-metallic nanoparticles (nanoalloys¹) exhibit size- and composition-dependent properties that can be exploited for various applications,² including catalysis,^3–6 plasmonics⁷ and chemical sensing.⁸ In catalysis, for example, the chemical activity is often attributed to a specific feature in the atomistic structure,^3–5 such as a particular facet³ or mixing/segregation pattern;⁴ and to guarantee the intended functionality the desired features must be thermodynamically stable or sufficiently long-lived in a given environment. This precondition naturally motivates the study of nanoalloy thermodynamics⁹ and calls for capabilities to manipulate finite-system analogues of a phase transition, which are usually associated with a relatively sudden loss of structural integrity around a particular temperature. In the present study we focus on solid–solid transitions, where one solidlike phase is supplanted by a distinctly different one prior to complete melting, and we explore how such transitions are affected when the atomic identity of a single substituent is varied. We consider various models and set three main objectives: (i) to qualitatively compare with previous experimental^10,11 and computational^12–19 studies examining impurity/dopant effects on cluster melting; (ii) to discern generic behaviour that could be exploited for tuning various thermal instabilities in nanoalloys; and (iii) to demonstrate how the harmonic superposition approach (HSA)^20,21 can be used to identify and explain impurity effects in atomic clusters at temperatures below melting.

By measuring the heat capacity as a function of temperature, Jarrold et al.^10,11 have shown that an aluminium impurity in size-selected Ga_N−1Al⁺ (N = 17, 19, 20, 30–33, 43, 46, 47) clusters has only a small effect on the melting behaviour,¹⁰ whereas a copper substituent in Al_N−1Cu⁻ (N = 49–62) generally causes a significant change.¹¹ From previous computational studies we also know that impurity effects on cluster melting depend on the impurity¹⁵ as well as the host cluster,¹⁸ suggesting that selective substitutional doping may be a feasible strategy for tuning the finite-temperature behaviour. Indeed, a single impurity can cause the melting temperature of a cluster to decrease¹³ in some cases and increase^14,16,17 in others. In particular, molecular dynamics simulations of Mottet et al.¹⁶ show the melting temperature of geometrically closed-shell Ag_N Mackay icosahedra²² increasing by up to 70 K when the central silver atom is substituted by nickel or copper, and this trend was correlated with impurity-dependent strain relaxation, which evidently favours the solidlike icosahedron over the liquidlike state.

The present study provides a broader evidence base for the behaviour described by Mottet et al.,¹⁶ demonstrating that it is not specific to melting. We show that a solid–solid transition can also be tuned via dopant-controlled stress/strain redistribution, provided: (i) the dopant-host atom-size mismatch is within a system-specific limit, which in turn can depend on the relative binding strength of the impurity; and (ii) one of the two competing phases is more inhomogeneously strained than the other, e.g. icosahedra competing with crystalline fcc fragments or Marks decahedra.²³ A weakly and sometimes moderately mismatched substituent will stabilise the more strained competitor, which will cause the transition temperature to shift in the appropriate direction by an amount that depends on the magnitude of atom-size mismatch and, to a lesser degree, the relative binding strength. In general, however, transmuting a single atom can also cause surprisingly intricate effects on a preexisting transition, often completely destroying it or producing multi-stage permutational isomerisation, which we discern for a variety of Lennard-Jones²⁴ (LJ) clusters and metal clusters modelled by the Gupta^25,26 potential.

All our thermodynamic analysis is based on the HSA,^20,21 which is a semi-quantitative approximation that has proved very useful in studies of atomic and molecular clusters.^21,26–28 The approximation entails coarse-graining the classical partition function into additive contributions from geometrically distinct local minima on the underlying potential energy landscape. The problem of thermodynamic sampling, which is notoriously difficult for solid–solid transitions,²⁹ is thus reduced to a one-off search for local minima. The search need not be constrained by any statistical distributions, because the approximate statistical weight of each minimum is given by a simple analytic form. The inherent separability and analyticity of the HSA make it particularly useful for developing a qualitative understanding of equilibrium structural transitions, and at sufficiently low temperatures, when anharmonic effects play a negligible role, the HSA is expected to be quantitatively accurate in the classical regime.

We revisit some general aspects of the HSA in section II, focusing on the origin and interpretation of heat capacity anomalies, and we consider illustrative examples based on generic two-state and three-state models. In section IIIA we define atomistic potentials to inform more complicated models for nanoalloys, while the procedure for obtaining a representative database of low-lying minima is outlined in section IIIB. Section IV contains all the results and discussion for selected Lennard-Jones (IVA) and Gupta (IVB) clusters, with the main conclusions summarised in section V. For completeness, in Appendix A we derive and carefully identify the inherent approximations of the HSA in the canonical ensemble.

II. General considerations

For a system of N classical atoms in volume V and at temperature T, the HSA amounts to approximating the canonical partition function as^20,21where m spans a representative database of M geometrically distinct local minima on the underlying potential energy landscape, U_m is the potential energy of the local minimum m, g_m is a degeneracy factor subsuming various entropic components, and k_B is the Boltzmann constant. Throughout this study m = 1 refers to the putative global minimum (i.e. the ground state). Appendix A provides a detailed derivation and discussion of the entropic factorwhere N_s is the number of s-type atoms (Σ_sN_s = N), κ = 3N − 6 is the number of vibrational degrees of freedom, h is the Planck constant, and o_m and [small nu, Greek, macron]_m are, respectively, the number of symmetry elements in the point group and the geometric mean normal mode vibrational frequency of m. Note that {U_m, [small nu, Greek, macron]_m, o_m}^M_m=1 are parameters whose values can be obtained from classical atomistic models (see section IIIA) or higher levels of theory.

Finite-system analogues of a phase transition can be identified from peaks (or other anomalies) in the heat capacity,

where the constant kinetic contribution has been omitted and the angle brackets 〈…〉_T indicate a canonical weighted average of the quantity within, i.e.withrepresenting the equilibrium occupation probability of each local minimum.

Within the limits of the harmonic approximation employed throughout this study, eqn (3) yields the vibrational component of the heat capacity. Note that C_V(T) = 0 when M = 1, but for M > 1 it is instructive to consider the heat capacity in a more explicit form:

where X = Δ₂₁/k_BT > 0, Δ_ij = U_i − U_j, ϕ_ij = Δ_ij/Δ₂₁, and ρ_ij ≡ g_i/g_j = ([small nu, Greek, macron]^k_jo_j)/([small nu, Greek, macron]^k_io_i) > 0. (Recall that the ground state has index 1.) It is now clear that the heat capacity depends on degeneracy ratios ρ_ij = ρ_ji⁻¹ = ρ_ikρ_kj and scaled energy differences, assuming j < i in the last equality. C^(M)_V tends to zero in the limits of T → 0 and T → ∞, so we expect at least one peak (sometimes referred to as a Schottky anomaly^30–32) in the range 0 < T < ∞ for finite M > 1.

For M = 2 the position X_max > 2 of the peak maximum is determined by the positive solution of

This equation is solved graphically in Fig. 1, illustrating how X_max > 0 depends on the degeneracy ratio ρ(= ρ₁₂), and also showing a second negative solution which is unphysical. Note that varying ρ or Δ₂₁ does not change the characteristic shape of the C_V peak.

Fig. 1 Graphical solution to eqn (6) and (9), demonstrating the effect of the degeneracy ratio ρ = g₁/g₂ in a two-state model.

The situation becomes more interesting for M = 3, because varying, ρ₂(= ρ₂₁), ρ₃(= ρ₃₁ = ρ₃₂ρ₂₁) and ϕ₃(= ϕ₃₁ = 1 + ϕ₃₂) affects the overall shape of the C_V curve and can yield two separate peaks, as illustrated in Fig. 2. The height and relative position of the peak(s) can be systematically adjusted, with the degeneracy ratios largely affecting the height, and the energy-level ordering/spacing primarily dictating the relative positions of the two maxima. Interestingly, bimodality seems to occur only when the higher energy gap is considerably wider than the low energy gap, and it is conceivable that an M-level spectrum with ever-increasing gaps between consecutive energy levels may have M − 1 peaks in the C_V curve. Although this limit is unlikely to be reached for atomistic systems, where high-energy minima are usually more densely spaced than low-energy minima, in section IVA we show that dopant-induced degeneracy splitting of the global minimum can result in more than one additional C_V peak at low temperatures.

Fig. 2 Heat capacity of a three-level system with different values of g₂, g₃ and ϕ₃, with Δ (i.e. the energy gap between the ground state and the first excited state) setting the energy scale. The default values are g₁ = g₂ = g₃ = 1 and ϕ₃ = 0.1. Note the logarithmic scale on the horizontal axis.

In section IV we consider databases containing M ≤ 10⁴ geometrically-distinct low-lying minima, with {g_m, U_m}^M_m=1 calculated from classical atomistic models. The overall shape of the resulting C_V plots is often qualitatively similar to that of a three-level system, but in some particular cases we encounter up to four identifiable peaks. Crucially, the HSA allows us to quantify the contribution of subsets of local minima to an overall C_V plot, which facilitates heuristic interpretation of various anomalies. A pronounced peak can often be associated with a transition from preferential occupation of a subset (α) of low-lying minima to a different subset (γ) of higher-lying minima.²⁶ We interpret the subsets α and γ as two different phases, with the number of distinct minima (i.e. |α| and |γ|) providing a measure of the landscape contribution to the configurational entropy³³ for each phase. Within this picture, the corresponding transition temperature T_s > 0 can be defined as the solution of Z_α(T) = Z_γ(T),²⁷ where represents a relative occupation probability. Equating these probabilities is equivalent to equating Helmholtz free energies, i.e. F_α(T) ≡ −k_BT ln Z_α, and determining T_s amounts to finding the solution(s) of

where the lowest-lying minimum in subsets α and γ has been factored out and

To show that a C_V peak does not always constitute a finite-system analogue of a phase transition, we again consider a two-state system with |α| = |γ| = ζ_α = ζ_γ = 1. Eqn (7) then yields a unique analytic solution:

ρ exp(X_s) = 1, (9)

where X_s = ΔU/(k_BT_s), and . The solution of (9) is compared to that of (6) in the inset of Fig. 1, illustrating that X_s < X_max and, hence, T_s > T_max for all ρ, with T_max denoting the temperature at the C_V maximum. While a C_V peak with T_max > 0 exists for all ρ, it constitutes a finite-system analogue of a phase transition with a physically meaningful T_s > 0 only when the high-energy phase has a higher entropy than the low-energy phase, i.e. ρ < 1. Also, T_max ≈ T_s > 0 for ρ ≪ 1, and in most of the atomistic models considered in section IVA we find the discrepancy between T_s and the C_V maximum is negligible.

The relatively simple expression in (9) is accurate only when the landscape component of the configurational entropy plays a negligible role. However, since we consider the contributions of both the landscape and the well vibrational entropies, we solve (7) numerically throughout this study (using the Powell hybrid method implemented in MINPACK³⁴). Furthermore, we cannot eliminate the possibility of (7) having multiple solutions, which would describe a reentrant phase with lowest free energy in multiple, disjoint regions of the temperature domain. However, none of the cases considered below exhibit any signs of reentrant behaviour.

III. Models and methods

A. Atomistic potentials

In this study we consider a model binary cluster containing atomic species A (the host) and B (the dopant), with each atom i identified by a label l_i ∈ {A, B}. The (binary) Lennard-Jones potential is then given bywhere r_ij is the distance between atoms i and j, ε_lilj = ε_ljli is the pair well depth, and σ_lilj = σ_ljli controls the equilibrium pair separation. In our analysis of LJ clusters (section IVA) we scale the length and energy by σ_AA and ε_AA, respectively, and treat the ratios σ_r = σ_AB/σ_AA and ε_r = ε_AB/ε_AA as adjustable parameters specifying the relative size and binding strength of the dopant atom. We also use the reduced temperature T* = k_BT/ε_AA.

For a more appropriate description of transition metal nanoalloys we employ the Gupta potential^25,26

where A_lilj, ξ_lilj, p_lilj, q_lilj and are system-specific parameters, and r_ij as in (10). We take the values fitted by Cleri and Rosato²⁶ for homonuclear (l_i = l_j) interactions, while the heteronuclear (l_i ≠ l_j) values are calculated using the Lorentz–Berthelot rules: by applying geometric averaging for A_lilj and ξ_lilj, i.e., and arithmetic averaging for the remaining parameters, i.e. p_lilj = (p_lili + p_ljlj)/2. These “average” parameters are not expected to be the most accurate values for specific nanoalloys, but they are appropriate for our present purpose, which is to explore qualitative trends with respect to the relative size and binding strength of the dopant atom, by analogy with the LJ model.

The LJ and Gupta potentials permit calculation of atomic level stresses³⁵ and, hence, the local (static) pressure P_i for each atom i,

where Ω_i is the effective atomic volume, F_ij is the force on atom i due to atom j, and r_ij = r_i − r_j. Note that Ω_i is often approximated by the Wigner–Seitz radius in the bulk ground state, irrespective of the particular cluster structure.^16,36

B. Databases of minima

Databases of local minima were harvested using generalised³⁷ basin-hopping^38,39 (GBH), with the details provided in ESI.† The database size is denoted by M, and the constituent local minima are indexed in the order of increasing potential energy. To rapidly survey in the neighbourhood of (ε_r, σ_r) = (1, 1) for doped LJ_N clusters, we took a simpler “enumerative” approach⁴⁰ analogous to that of Cao et al.¹¹ Starting from a database for a given homoatomic LJ_N system, split into sets α_N and γ_N, we systematically generated α′_N(ε_r, σ_r) and γ′_N(ε_r, σ_r) by enumerating all the N permutations for each structure (referred to as the “parent”) in α_N and γ_N, respectively, with one atom labelled as the dopant. The resulting configurations in α′_N(ε_r, σ_r) and γ′_N(ε_r, σ_r) were then requenched, all duplicates were removed, and the two sets of minima were fed into eqn (7) to calculate . While this approach avoids the cost of running GBH for each (ε_r, σ_r) independently, it may fail to locate new low-energy minima that can emerge when the mismatch in atom size or binding strength exceeds a certain (unknown) threshold. Hence, we restricted ourselves to a relatively small window in the two-parameter space, i.e. 0.97 ≤ σ_r ≤ 1.03 and 0.65 ≤ ε_r ≤ 1.5, and we checked the validity of α′_N(ε_r, σ_r) and γ′_N(ε_r, σ_r) at the four corners of this window by comparing with an independent GBH run.

IV. Results and discussion

A. Lennard-Jones clusters

We start by reviewing known^26,27,41,42 solid–solid transitions in LJ_N clusters with N = 31, 38, 75 and 98. Fig. 3 shows how these transitions are manifested in a heat capacity plot: either as an isolated peak (LJ₃₁, LJ₇₅ and LJ₉₈) or as a shoulder of a higher temperature peak (LJ₃₈). Plotting the occupation probabilities reveals that the low-temperature C_V anomaly in each case correlates with destabilisation of the global minimum (GM), whose occupancy drops below 0.5 at a system-specific reduced temperature and continues to decrease. To determine for each N, we solve (7) with the subset α containing just the corresponding GM, and the subset γ containing a particular number of consecutive higher-lying minima. The subset size |γ| is determined by first including just the second-lowest minimum in γ, then adding higher-lying minima in the database (of size M) one by one, in the order of increasing potential energy, until the numerical solution to eqn (7) no longer changes. The smallest |γ| required for acceptable convergence and the corresponding values are listed in Table 1 for all N considered, and we checked that these temperature values change by less than 1% if eqn (23) is used instead of (2).

Fig. 3 Occupation probabilities (P_occ) of the ten lowest-energy minima and heat capacities (C_V) plotted versus reduced temperature for (a) LJ₃₁, (b) LJ₃₈, (c) LJ₇₅ and (d) LJ₉₈. The occupation probabilities are stacked on top of each other in the order of increasing potential energy and colour-coded (see the key). The two lowest-energy minima in each case are shown for illustration. Converged C_V plots are represented by a solid black curve, while the effect of excluding the global minimum is shown by a dashed curve. Convergence of the solid–solid peak is demonstrated by (thinner) intermediate C^(M)_V(T*) curves, with the corresponding size (M) of the minima database indicated by an integer. Vertical dashed lines mark the transition temperature obtained by solving eqn (7) numerically, and the horizontal dashed lines correspond to P_occ = 0.5.

Table 1 Transition temperature for selected LJ_N clusters (calculated using 1 + |γ| lowest-lying minima) and the range of local pressure values (multiplied by a volume) for the lowest-lying minimum in subsets α (with |α| = 1) and γ

N		|γ|	P ^(α)min	P ^(α)max	P ^(γ)min	P ^(γ)max
10	0.0271	3	−2.693	13.875	−4.491	18.210
38	0.1209	16	−0.469	1.248	−1.830	13.461
75	0.0824	17	−0.734	2.564	−3.413	29.501
98	0.0035	1	−1.519	2.904	−3.446	30.185

---

The GM of LJ₃₁ is an incomplete Mackay icosahedron,²² and the structural transition in Fig. 3a can be described as surface polymorphism: the incomplete outer shell gains access to multiple configurations, one of which is anti-Mackay (the second-lowest minimum). For LJ₃₈, the dominant C_V peak in Fig. 3b has also been interpreted as a Mackay-to-anti-Mackay transition,⁴¹ while the low-temperature shoulder corresponds to the octahedral GM being supplanted by several incomplete icosahedra. Note that it is predominantly the third- and fifth-lowest minima that supplant the LJ₃₈ GM, with the second- and fourth-lowest minima playing a negligible role (due to relatively low vibrational entropy). This observation shows that the ordering by potential energy may not necessarily reflect the precedence/occupancy at finite temperatures, which is perhaps not surprising since it is the free energy that becomes the determining factor. More importantly, it shows that one should not assume that local minima lying immediately above the GM will be statistically important.

For LJ₇₅, the C_V peak in Fig. 3c corresponds to the Marks decahedral²³ GM giving way to numerous incomplete Mackay icosahedra.²² Although this particular transition has been discussed in some detail,^27,43 we point out an apparent similarity with structural transitions occurring in certain metal clusters during solid–liquid phase coexistence.⁴⁴ In the latter systems the formation of a significant liquidlike fraction at constant energy coincides with a sudden structural transition in the solidlike core. The present analysis of LJ₇₅ suggests that the onset of surface polymorphism at constant temperature can also coincide with abrupt structural transformation of an entire cluster.

The solid–solid transformations in LJ₃₁, LJ₃₈ and LJ₇₅ can be classified as “one-to-many”,²⁸ because the GM is supplanted by multiple isomers with greater vibrational entropy. In contrast, Fig. 3d reveals a “one-to-one” transformation for LJ₉₈: the two lowest-lying minima are sufficient to converge the low-temperature C_V peak, while the occupation probabilities of higher-lying minima remains negligible until significantly higher temperatures. The same two-state model has been shown to reproduce the low-temperature heat capacity of LJ₉₈ calculated using replica-exchange Monte Carlo quite well.⁴² Note that the transition temperature obtained using (7) is noticeably above the peak maximum, and the origin of this discrepancy was discussed in section II (Fig. 1 inset).

We now survey as a function of ε_r and σ_r using the “enumerative” procedure described in section IIIB. The results (interpolated using a cubic spline) are plotted in Fig. 4, illustrating how is affected by small changes in the dopant characteristics. Although the accuracy of these “heat-maps” is expected to deteriorate when moving away from (ε_r, σ_r) = (1, 1) in the two-parameter space, the proximity of dark blue regions (where ) to the point (ε_r, σ_r) = (1, 1) indicates how sensitive (or robust) a given transition is to the dopant-host mismatch. We see that is generally more sensitive to size mismatch (σ_r) than relative binding strength (ε_r), which seems intuitive and consistent with the only available experimental data on the melting of doped metal clusters.^10,11 What is less intuitive is that the transitions in LJ₃₈ and LJ₇₅ are significantly more sensitive to a smaller dopant (σ_r < 1), whereas LJ₃₁ is more sensitive to a larger dopant (σ_r > 1). We address this particular issue below. Another predictable observation is that generally increases with ε_r > 1, and perhaps it eventually catches up with the system-specific melting temperature, but we shall not focus on this scenario. It is also worth mentioning that the “heat-map” for LJ₉₈ (not shown) resembles those of LJ₃₈ and LJ₇₅, exhibiting a similar pattern but within a significantly narrower region around (ε_r, σ_r) = (1, 1).

Fig. 4 Variation of the solid-solid transition temperature with the relative size (σ_r) and binding strength (ε_r) of the dopant atom for 31 (top), 38 (middle) and 75-atom (bottom) clusters.

It is instructive to consider the effect of varying σ_r while keeping the binding strength fixed at ε_r = 1,⁴⁰ in which case increasing or decreasing the dopant size lowers (see Fig. 5). We link this behaviour to the homoatomic GM exhibiting a more homogeneous distribution of local pressure than the higher lying icosahedral minima (see Table 1). The centre of an icosahedral motif is under a particularly high compressive stress, which can be relieved by substituting the central atom with a smaller one.^16,40 Hence, of doped LJ₃₈ and LJ₇₅ rapidly drops to zero in Fig. 4 when σ_r is below a critical value σ^†_r(ε_r) < 1, which in turn can be regulated by varying ε_r.

Fig. 5 Scaled versus dopant relative size (σ_r) for selected LJ_N clusters. Dashed and dotted curves represent data with the vibrational and landscape components of the configurational entropy, respectively, set to the values for the homoatomic cluster (σ_r = 1).

Fig. 5 clearly shows that in doped LJ₃₈, LJ₇₅, and LJ₉₈ the transition temperature is more sensitive to σ_r < 1 than it is to σ_r > 1,⁴⁰ and in the latter case the dopant energetically favours a surface site. This asymmetric sensitivity can also be explained by considering the local pressure: the difference in maximal expansive (negative) pressure between the competing morphologies is significantly smaller in magnitude than the difference in maximal compressive (positive) pressure, so a larger dopant has a smaller effect on . The situation is in some sense reversed for LJ₃₁, where is more sensitive to a larger dopant. In this case the competing morphologies share the same icosahedral core and are therefore similarly affected by a smaller dopant. However, a larger dopant atom occupying a surface site evidently has a greater destabilising effect on the Mackay capping layer, which seems intuitive since the anti-Mackay capping is less densely packed (see snapshots in Fig. 3a) and thus better able to accommodate expansive strain.

Rationalising the trends in Fig. 5 (and 4) solely in terms of local pressure is justified only if the entropic component of the impurity effect can be neglected. To check the role of vibrational entropy, we considered the case where minima in α′(σ_r) and γ′(σ_r) have [small nu, Greek, macron] set to the value for their “parent” (recall section IIIB) in α and γ, respectively. This artificial shift undoes the changes in vibrational entropy induced by adjustments in σ_r, and we find that it does not affect the qualitative nature of the trends in Fig. 5. Effects of the landscape contribution to configurational entropy can also be (largely) eliminated by restricting α′(σ_r) and γ′(σ_r) to include only biminima^37,45—configurations whose energy cannot be improved by interchanging any two atoms (and then requenching)—with the corresponding point group order manually changed to that of the “parents” so as to preserve continuity at σ_r = 1. This manipulation has a more drastic effect on the curve, as shown for LJ₃₈ and LJ₇₅ in Fig. 5, but the overall qualitative picture remains unchanged: is still highest at σ_r = 1, and the mismatch-induced lowering of is still more gradual for σ_r > 1. Hence, the effect of σ_r on is determined primarily by changes in the potential energy, at least for the LJ clusters considered here, with entropy playing only a secondary (though quantitatively important) role. This conclusion suggests that impurity atoms energetically favouring interior sites within a cluster will generally have a stronger effect on structural transitions, because a higher coordinated impurity will have a greater share of the binding energy. However, solid–solid transitions that do not involve core restructuring may not follow this rule, as evidenced by LJ₃₁.

Fig. 4 and 5 illustrate systematic shifts in for a geometric solid–solid transition, and Fig. 6a shows the shift for LJ₇₅ happening in unison with the computed using (7). (Similar consistency has also been demonstrated⁴⁰ for LJ₃₁ and LJ₃₈.) Furthermore, Fig. 6a exhibits additional low-temperature peaks that also change with σ_r. Analogous C_V features have been reported previously¹⁹ and attributed to dopant delocalisation prior to melting. Here we also interpret these features as indicators of non-degenerate permutational isomerisation (see below), and the fact that they change systematically with σ_r indicates potentially tunable behaviour. Interestingly, Fig. 6b shows that eliminating σ_r-induced variation in the configurational entropy accentuates the main C_V peak and preserves the direction in which it shifts in response to σ_r, while the features associated with dopant delocalisation vanish completely. It would be interesting to investigate if the result will be the same for multiple impurity atoms, or whether some C_V anomalies due to permutational isomerisation will be retained, but we will not pursue this issue further here.

Fig. 6 Heat capacity of doped LJ₇₅ clusters for different values of σ_r, calculated using either (a) all local minima or (b) just biminima. Dashed vertical lines mark the value of calculated using (7).

Fig. 7 focuses on the precursor anomalies indicated in Fig. 6, revealing two separate peaks for doped LJ₇₅ with σ_r = 0.992. This bimodality is apparently due to the dopant delocalising in two separate stages. The first stage corresponds to partial delocalisation of the dopant from the central site to the two adjacent ones along the five-fold symmetry axis. This interpretation is based on the fact that the lowest-temperature C_V peak in Fig. 7 is reproduced by the two lowest-lying homotops with D_5h and C_5v point group symmetry. Note that eqn (7) yields for this two-state subsystem, while the C_V peak maximum is more in line with the first inflection point (near T* = 0.001) in the GM occupation probability. The disparity between the C_V maximum and is consistent with the degeneracy ratio (recall Fig. 1 inset and section II) being close to unity, i.e. ρ = ([small nu, Greek, macron]^k₂o₂)/([small nu, Greek, macron]^k₁o₂) ≈ 0.5 in this case. The second precursor peak in Fig. 7 converges when all ten distinct permutations are considered, and the peak maximum at around T* = 0.01 is in better agreement with the solution of eqn (7), i.e. if the transition is between sets of two and eight minima. Hence, the second delocalisation stage here corresponds to the dopant gaining access to all the remaining sites within the underlying framework.

Fig. 7 Occupation probabilities, stacked and colour-coded as in Fig. 3, and (partial) heat capacity plots focusing on the low-temperature precursor anomalies in Fig. 6 for σ_r = 0.992. The five lowest-energy minima are shown for illustration, viewed down the principal axis (top row) and along the orthogonal plane (bottom), with the larger red sphere representing the dopant.

The origin of two-step permutational isomerisation in doped LJ₇₅ is analogous to the generic three-state model discussed in section II. Fig. 8 shows that atom-size mismatch (σ_r ≠ 1) splits the degeneracy of the homoatomic (σ_r = 1) energy levels, and for σ_r < 1 the degeneracy of the homoatomic GM is split into two well-defined bunches. A narrow energy gap between the two energy levels in the low-lying bunch (i.e. the first two isomers illustrated in Fig. 7) and the significantly wider gap between the two bunches give rise to two C_V peaks at low temperature. In contrast, the homoatomic GM degeneracy splitting for σ_r > 1 yields a more uniformly dispersed subspectrum, without anomalously large gaps or bunching, producing a very broad (and hardly noticeable) C_V peak.

Fig. 8 The potential energy spectrum of low-lying minima (top) and low-temperature heat capacity map (bottom) for doped LJ₇₅ as a function of the dopant relative size σ_r. In the top diagram, energy levels emanating from the same degenerate homoatomic (σ_r = 1) minima are coloured the same. The bottom diagram highlights the low-temperature (T* < 0.25) C_V peak(s) corresponding to permutational isomerisation.

The main solid–solid peak for each σ_r in Fig. 6 is due to the wide gap at the low end of the energy spectrum for σ_r = 1, and Fig. 8 shows that introducing dopant-host mismatch (i.e. increasing |1 − σ_r|) effectively narrows this gap and eventually closes it completely. This gap closing has the effect of smearing and eventually eliminating the C_V peak associated with the transformation from decahedral to icosahedral geometry. This observation reinforces the point made by Bixon and Jortner,⁴⁶ that pronounced gaps in the energy spectrum play a central role in determining finite-system analogues of a phase transition, while a gapless spectrum with uniformly distributed energy levels will not produce sharp peaks in the heat capacity. We expect that multiple impurities will result in more intricate degeneracy splittings than in Fig. 8 and, hence, will likely produce a more uniform energy spectrum and a more smeared C_V curve. This expected smearing would be analogous to the broadening of the specific heat divergence at a critical temperature,⁴⁷ and it should be distinguished from the broadening arising from finite-size effects.^48,49

To conclude our discussion of LJ_N clusters, we reiterate that geometric solid–solid transformations for N = 38, 75 and 98 are all affected very similarly when a substituent atom is transmuted. This similarity stems from the fact that these transformations correspond to a highly symmetric non-icosahedral GM being supplanted by numerous incomplete icosahedra. In the following section we consider model nanoalloys where the icosahedral order is less prevalent (for known reasons^50–53), and the diversity of solid–solid transitions is considerably richer. These examples correspond to the discrete values provided for us by the periodic table.

B. Gupta clusters

We used the HSA to explore solid–solid transitions in Ag_N, Al_N, Au_N, Cu_N, Ni_N, Pd_N, and Pt_N clusters with N = 38 and 75, and in all these cases the GM has the same atomic structure as the LJ_N counterpart, but the structure and energetic ordering of higher lying minima varies among metals. We found cases with multimodal C_V plots, with the structural interpretation of each peak also varying among metals, but we leave the details for a separate study.

A solid–solid transition most closely resembling its LJ counterpart was found in Ni₇₅, where the Marks decahedral²³ GM is supplanted by numerous incomplete Mackay icosahedra²² at k_BT ≈ 31 meV (see Fig. 9). The insight gained from LJ_N clusters leads us to expect that substitutional doping will lower the transition temperature of Ni₇₅, provided the dopant binding is weaker than Ni–Ni binding. A database of around 10⁴ low-lying minima (harvested using the GBH procedure described in ESI†) for Ni₇₄X clusters (with X = Au, Ag, Al, Cu, Ni, Pd, Pt, Pb) largely confirms this expectation. Fig. 9 shows that the most closely size-matched dopant (Cu) yields the smallest downward shift, while the most size-mismatched dopant (Pb) yields the largest shift. Note that here we define the overall dopant-host mismatch based on the relative position of different atomic species in Fig. 10, which illustrates where various metals lie in an empirical two-parameter space and can provide some guidance when selecting a dopant for a given host. Also note that the shift in Fig. 9 is less clear-cut for substituents with intermediate atom-size mismatch, indicating that additional factors may also be involved here.

Fig. 9 Top: Occupation probabilities (stacked on top of each other) versus temperature for Ni₇₅, based on the 40 lowest-energy minima, with 14 of them represented individually by a color (see key), and the remaining ones lumped into the “residual” (black). The horizontal dashed line marks P_occ = 0.5, while the vertical dashed line indicates the temperature at which the occupation probability of the GM crosses this value. Bottom: Heat capacity plots for Ni₇₄X clusters, showing how the solid–solid peak tends left as the dopant (X) atom size increases. The dashed black curve corresponds to the homoatomic Ni₇₅ with the decahedral GM removed from the minima database.

Fig. 10 Scatter plot of cohesive energy E_coh (per atom) and effective atomic radius R_eff for various 12-coordinated (i.e. fcc or hcp) metals.^26,54

The Pb-induced downward shift of the C_V peak by k_BΔT ≈ 15 meV (ΔT ≈ 170 K) is reinforced by the considerably smaller cohesive energy of Pb in comparison to Ni. However, the Pt-induced downward shift indicates that atom-size mismatch is the more important factor, because the considerably higher cohesive energy of Pt (relative to Ni) would push the C_V peak up were it not for the larger effective size of Pt. The shift also correlates with local atomic-level pressure, by analogy with LJ₇₅: a larger-sized dopant stabilises the incomplete icosahedral framework, because it contains sites with higher maximal expansive stress. The lowest value of the quantity defined in (12) is −0.55 eV for the Marks decahedral GM and −1.10 eV for the lowest-lying icosahedron-based isomer of Ni₇₅, and in both structures sites with the maximal expansive pressure are localised at (or near) pentagonal vertices.

The C_V plots in Fig. 9 also show broad peaks due to permutational isomerisation in Ni₇₄Ag and Ni₇₄Pt. These precursor peaks arise just from the two lowest-lying minima in each case, which are illustrated in Fig. 11. Although the situation is analogous to that in doped LJ₇₅ with σ_r > 1 (recall Fig. 8), it is interesting that the two competing permutational isomers are different for each dopant, which can be explained by the disparity in the cohesive energies of Ag and Pt. This observation shows that subtly different permutational-isomerisation transitions can occur in clusters with the same geometric framework.

Fig. 11 The two lowest-energy minima for (a) Ni₇₄Ag and (b) Ni₇₄Pt Marks decahedra,²³ viewed down the principal axis (top) and along the orthogonal plane (bottom). The larger red sphere represents the dopant, and the host atom size is varied for clarity.

Metal clusters modelled by a many-body potential generally exhibit a diminished propensity for icosahedral order^50–53 compared to Lennard-Jones clusters, with Au₅₅ providing the most striking example.^51–53 Although this size corresponds to a “magic” number (i.e. geometrically closed-shell) Mackay icosahedon,²² the GM is actually a defective fcc/hcp structure with C₁ point group symmetry. Fig. 12 shows its occupation probability diminishing relative to several amorphous structures with higher vibrational entropy, with the transition marked by a pronounced C_V peak at k_BT ≈ 8 meV. Although the six lowest-energy minima are all predominantly fcc-like and therefore may be regarded as constituents of one phase, the combined occupation probability of isomers 2 to 6 never exceeds 0.07. Isomers 7, 8, 9 and 14 collectively become the most populated at k_BT_s ≳ 8 meV, constituting a phase that can be characterised (using common-neighbour analysis^28,49,55) as amorphous with partial icosahedral order.

Fig. 12 Occupation probabilities (stacked on top of each other) and heat capacity plots for Au₅₅, with three of the most populated structures illustrated.

We note that Garzón et al.^52,53 have used the Rosato, Guillope and Legrand (RGL)⁵⁶ parametrization of the Gupta potential and suggested an amorphous GM structure for Au₅₅. However, we find that reoptimising the fcc GM illustrated in Fig. 12 with the reduced RGL parameters (i.e. with ξ set to unity as in ref. 53) yields −52.65827, which is slightly lower than that of the energy of the lowest-lying amorphous isomer (−52.65684). It is also worth noting that in both cases the homoatomic Mackay icosahedron²² is not even in the 10⁴ lowest-lying minima for Au₅₅.

The competing structures contributing to the C_V peak in Fig. 12 exhibit a similar homogeneous distribution of local pressure, so the stress-redistribution argument employed in section IVA and for Ni₇₅ does not apply in this case, and we cannot say anything about the potential tunability. In fact, we find that the low-temperature peak essentially vanishes when one Au atom is substituted by a closely-matched one (e.g. Ag or Al shown in Fig. 13b), because it energetically destabilises fcc order by promoting one of the amorphous isomers to the GM. Similar reordering occurs for Pb and Pt dopants, and the pronounced low-temperature peak in each case indicates a transition of a distinctly different kind: one where an amorphous GM gives way to other amorphous structures with higher vibrational entropy, whilst the low-lying fcc/hcp isomers are statistically insignificant.

Fig. 13 C _V plots for Au₅₄X clusters, with the homoatomic Au₅₅ represented by a dashed black curve. Cluster snapshots illustrate the pristine I_h structure (left), a single rosette-like defect (middle), and a double paired rosette defect (right). The dopant is occupying the central site (not visible), while atoms belonging to the hexagonal rosette rings are depicted in red.

The interpretation of the low-temperature peak again changes when one Au atom is substituted with Pd, Cu or Ni, which all have a smaller effective atom size than Au. Each one of these impurities promotes the Mackay icosahedron²² to the GM, and the low-temperature C_V peak corresponds to partial amorphization of the icosahedron surface via the rosette-like mechanism,⁵⁷ with a fivefold vertex (or two adjacent ones) transforming into a hexagonal ring, while the impurity remains at the central site. This type of solid-solid transition occurs because surface amorphised icosahedra have more vibrational entropy. It is interesting that decreasing the dopant atom size (i.e. Pd → Cu → Ni) energetically stabilizes the ideal icosahedron and just the two partially amorphised variants illustrated in Fig. 13a, with the perfect icosahedron stabilised relatively more than the amorphised variant(s), hence the peak in Fig. 13a shifts to the right by as much as k_BΔT ≈ 20 meV (ΔT ≈ 230 K). This trend demonstrates that selective substitutional doping can help to suppress surface amorphisation, which could be particularly important for catalysis, and it is conceivable that similar trends may be found for other types of surface reconstruction/roughening in nanoalloys. The behaviour in X-doped Au₅₄X clusters also leads to two somewhat counterintuitive points: (i) a cluster surface transition can be tuned by substituting a central atom; and (ii) a highly mismatched impurity can lead to more predictable behaviour than a closely-matched one.

V. Summary and conclusions

We have demonstrated how the harmonic superposition approach^20,21 can be used to explain the effects of substitutional doping on solid–solid transitions in atomic clusters. We first considered a selection of Lennard-Jones (LJ_N<100) clusters with a preexisting transition and analysed the effects of small changes in the relative atom size and binding strength of the impurity. We then focused on Ni₇₄X and Au₅₄X clusters (with X = Au, Ag, Al, Cu, Ni, Pd, Pt, Pb), modelled by a many-body Gupta potential, to further illustrate that a single impurity atom can produce systematic trends and/or surprisingly intricate (yet predictable) effects on a solid–solid transition in nanoalloys. Our analysis corroborates earlier studies of impurity effects^16,17 and sheds new light on the tunability of morphological transitions. We now summarise our main conclusions, which we believe may be useful for rational design of nanoalloy catalysts.

A solid–solid transition temperature, when defined by equating the free energies of two competing phases, is generally more sensitive to the relative size than to the binding strength of the impurity atom. This result is consistent with previous experiments^10,11 and calculations.^{12–16,18,19} More importantly, we have shown that a given solid–solid transition supports only a limited range of atom-size mismatch, and exceeding the system-specific limits will destroy the transition. Hence, a transition temperature can be tuned only by impurities with mismatch magnitude below a certain tolerance, which is generally larger for surface (as opposed to interior) impurities. The tolerance does not exhibit a simple dependence on system size, but rather is determined by details of the underlying potential energy landscape.

From the energy landscape perspective, small changes in the dopant characteristics cause an appreciable and systematic redistribution of local minima along the energy axis, while the effect on the corresponding local curvatures is relatively insignificant. Hence, impurity effects on solid–solid transitions are largely manifested in the discrete energy spectrum of low-lying minima (i.e. metastable states), which describes the main energetic factors as well as the landscape contribution to configurational entropy.³³

In certain tunable cases we found that dopant effects can also be linked to specific features in the distribution of atomic-level stresses. This correlation was discussed in the context of melting by Mottet et al.,¹⁶ who showed that the solid–liquid transition temperature of Ag clusters can be raised by a smaller impurity. Our analysis extends these results by showing that both larger and smaller impurities can lower a solid–solid transition temperature. These seemingly opposing trends are consistent with one general principle: a size mismatched impurity will generally favour a more inhomogeneously stressed phase. This principle can help predict the direction in which a transition temperature will shift after doping, provided the atomic-level stress in structures representing the two competing phases are known.

Interestingly, we found that a single impurity atom can give rise to more than one additional anomaly in the low-temperature heat capacity, indicating that dopant delocalisation within a particular geometric framework can proceed via two separate stages (and possibly more in larger systems). This result constitutes a simple analogue of the behaviour reported by Rubinovich et al.³² for model Pd–Cu nanoalloys, where the solid-solution-like transition occurred in multiple stages of partial mixing/segregation staggered over a finite temperature range.

From a modelling perspective, the difficulty of examining multiple-impurity effects will grow rapidly with the number of impurities, because the set of distinct local minima required for the harmonic superposition approach will be significantly greater than those used here. However, we believe this difficulty can be (at least partially) overcome by systematically reducing the database of minima to a smaller representative set containing only biminima.⁴⁵ In the present study we have shown that biminima alone can retain the necessary features that give rise to geometric structural transitions, and it would be interesting to examine to what extent permutational transitions are preserved by biminima when there is more than one impurity atom.

Appendix A: HSA derivation

Consider a system of N potentially distinguishable particles, each described by Cartesian coordinates r ∈ [Doublestruck R]³ (restricted to some finite volume V) and velocities ṙ ∈ [Doublestruck R]³. The canonical partition function Z(N, V, T) can be written aswhere Ω represents the region in 3N-dimensional configuration space imposed by the fixed volume V, x = (r₁, …, r_N), ẋ = (ṙ₁, …, ṙ_N), T is the temperature,is the classical Hamiltonian, U(x) is the potential and M is a diagonal matrix containing the atomic mass for each kinetic degree of freedom. Note the non-standard integration with respect to velocities ẋ, rather than momenta p = Mẋ, hence the Jacobian |M⁻¹|. Also note that the range of integration for each velocity component is (−∞, + ∞).

Now consider a local approximation to U(x) based on the Taylor expansion about a local minimum x_m:

where U_m = U(x_m); the first derivatives vanish; and the superscript t indicates the matrix transpose. The Hessian matrix H^(m) can be diagonalised by switching to normal mode coordinates using the following transformation:where A_ij is a 3N × 3N orthogonal matrix satisfying [A^tA]_ij = δ_ij andwith λ_i the i'th eigenvalue of the mass-weighted Hessian , and δ_ij the Kronecker delta. The eigenvector associated with λ_i is given by the i'th column of A. Also, note that it is not necessary to transform the velocities ẋ in unison with the coordinates x, but doing so is useful for separating out the kinetic energy contributions associated with particular normal modes. Now, substituting Ũ_m(x) for U(x) in (14) and transforming the coordinates according to (15) yieldswhere the subscript m identifies the local minimum about which the approximation is constructed. Substituting (17) for in (13) and transforming the differentials, i.e., yieldswhere i spans κ vibrational modes with λ^(m)_i > 0, and j spans the residual 3N − κ modes with λ^(m)_j = 0, if they exist. Note that the potential energy of an isolated cluster or a molecule with N > 2 will yield six eigenvalues that are zero at a stationary point: three due to translational invariance and the other three due to rotational invariance. However, all the degrees of freedom can become vibrational, e.g. in the presence of an external field, in which case κ = 3N and q^res_m = 1.

Now, each configuration integral (with respect to a component of y) ought to be evaluated over a finite range l^(m)_i/j for the catchment basin of m, which cannot be easily determined. For vibrational modes, it is conventional to extend the range of integration over all space, i.e. l^(m)_i → (−∞, +∞); and, using the standard result for a > 0, it becomes possible to evaluate the κ double integrals analytically to obtain

which represents the local vibrational entropy withthe geometric mean vibrational frequency for the local minimum m. The 3N − κ integrals defining q^res_m in (18) are usually separated into rotational and translational components, i.e. q^res_m = q^rot_mq^trans_m, and the two contributions are usually taken as

q^trans_m = Vh⁻³M^3/2(2πk_BT)^3/2, (21)

q^rot_m = 8π²h⁻³|I_m|^1/2(2πk_BT)^3/2, (22)

where M is the sum of atomic masses and |I_m| is the determinant of the inertia tensor. Note that the expression for q^trans_m would be exact if periodic boundary conditions were applied, but for a finite container it is actually an approximation, because the centre of mass of an extended body will be confined to a volume smaller than V. Also, for local minima with point group symmetry, q^rot_m will include contributions from superimposable configurations that are related by a proper rotation in the group.^21,58

A global approximation to the partition function (13) is constructed by superposing the local approximations for geometrically-distinct local minima, yielding eqn (1) with the degeneracy factor

g_m = n_mq^trans_mq^rot_mq^vib_m, (23)

where n_m accounts for non-superimposable symmetry-equivalent versions of m and is given by^21,58with s spanning different atomic species (satisfying N = Σ_sN_s), the factor of 2 explicitly accounting for inversion isomers, and o_m specifying the point group order. Division by o_m corrects for rotational and inversion isomers that are superimposable, as they also fall into the set of Π_sN_s! permutational isomers.

Note that the temperature dependence in q^vib_m, q^rot_m and q^trans_m has no effect on thermodynamic quantities due to cancellations when computing ensemble averages using eqn (4). Further, recent studies^28,59 show that the moment-of-inertia terms (i.e. |I_m|) in q^rot_m also largely cancel out, which is why we use the simpler expression for the effective degeneracy stated in (2).

Finally, we restate the three main sources of error intrinsic to the HSA, which have usually been discussed in the microcanonical (as opposed to canonical) context.⁶⁰ Firstly, errors will arise when the database of local minima is not fully representative, but these errors can in principle be eliminated by systematically extending the database. The second source of error is due to the neglect of anharmonic terms in the local approximations of the potential, which can also be systematically corrected by incorporating higher-order terms in the Taylor series expansion. Finally, errors can arise from the overlap between different local approximations in the configuration domain. This particular error is introduced when the finite limits of integration in (18) are extended over an infinite range, which is, in fact, necessary for analytic evaluation. In principle, each configuration integral should be evaluated over the basin of attraction for the local minimum, but this cannot be done analytically. Hence, correcting for the overlap error would inevitably require sacrificing much of the analytic utility of the HSA.

Acknowledgements

This work was funded by the ERC and EPSRC grant EP/J010847/1. BEH also acknowledges the Gates Cambridge Trust for financial support. The authors are very grateful to Dr Aleks Reinhardt for a critical reading of the manuscript, and DS would like to thank Marco Eckhoff with Dr Andy Ballard for fruitful discussions of the harmonic superposition approach.

References

1. R. Ferrando, J. Jellinek and R. L. Johnston, Chem. Rev., 2008, 108, 845.
2. Nanoalloys: from fundamentals to emergent applications, ed. F. Calvo, Elsevier, 2013.
3. K. Zhou and Y. Li, Angew. Chem., Int. Ed., 2012, 51, 602.
4. D. Wang, H. L. Xin, R. Hovden, H. Wang, Y. Yu, D. A. Muller, F. J. DiSalvo and H. D. Abruña, Nat. Mater., 2013, 12, 81.
5. S. Shan, V. Petkov, B. Prasai, J. Wu, P. Joseph, Z. Skeete, E. Kim, D. Mott, O. Malis, J. Luo and C.-J. Zhong, Nanoscale, 2015, 7, 18936.
6. W. Luo, M. Sankar, A. M. Beale, Q. He, C. J. Kiely, P. C. Bruijnincx and B. M. Weckhuysen, Nat. Commun., 2015, 6, 6540.
7. V. Amendola, S. Scaramuzza, S. Agnoli, S. Polizzi and M. Meneghetti, Nanoscale, 2014, 6, 1423.
8. S. Tong, Y. Xu, Z. Zhang and W. Song, J. Phys. Chem. C, 2010, 114, 20925; X. Li, J. Yao, F. Liu, H. He, M. Zhou, N. Mao, P. Xiao and Y. Zhang, Sens. Actuators, B, 2013, 181, 501; C. Wadell, F. A. A. Nugroho, E. Lidstrm, B. Iandolo, J. B. Wagner and C. Langhammer, Nano Lett., 2015, 15, 3563.
9. F. Calvo, Phys. Chem. Chem. Phys., 2015, 17, 27922.
10. C. M. Neal, A. K. Starace and M. F. Jarrold, J. Phys. Chem. A, 2007, 111, 8056.
11. B. Cao, A. K. Starace, C. M. Neal, M. F. Jarrold, S. Núñez, J. M. López and A. Aguado, J. Chem. Phys., 2008, 129, 124709.
12. J. C. Shelley, R. J. L. Roy and F. G. Amar, Chem. Phys. Lett., 1988, 152, 14; I. Garźon, X. Long, R. Kawai and J. Weare, Chem. Phys. Lett., 1989, 158, 525; I. L. Garźon, X. P. Long, R. Kawai and J. H. Weare, Z. Phys. D: At., Mol. Clusters, 1989, 12, 81.
13. C. Hock, S. Straßburg, H. Haberland, B. von Issendorff, A. Aguado and M. Schmidt, Phys. Rev. Lett., 2008, 101, 023401; A. Lyalin, A. Hussien, A. V. Solov'yov and W. Greiner, Phys. Rev. B: Condens. Matter, 2009, 79, 165403.
14. U. Ojha, K. G. Steenbergen and N. Gaston, J. Chem. Phys., 2013, 139, 094309.
15. A. Aguado, L. E. González and J. M. López, J. Phys. Chem. B, 2004, 108, 11722.
16. C. Mottet, G. Rossi, F. Baletto and R. Ferrando, Phys. Rev. Lett., 2005, 95, 035501.
17. D. Cheng, S. Huang and W. Wang, Phys. Rev. B: Condens. Matter, 2006, 74, 064117.
18. P. Chandrachud, K. Joshi and D. G. Kanhere, Phys. Rev. B: Condens. Matter, 2007, 76, 235423.
19. N. Quesada and G. E. Moyano, Phys. Rev. B: Condens. Matter, 2010, 82, 054104.
20. F. H. Stillinger and T. A. Weber, Phys. Rev. A, 1982, 25, 978; D. J. Wales, Mol. Phys., 1993, 78, 151; G. Franke, E. R. Hilf and P. Borrmann, J. Chem. Phys., 1993, 98, 3496; F. Calvo, J. P. K. Doye and D. J. Wales, J. Chem. Phys., 2001, 115, 9627.
21. D. J. Wales, Energy landscapes, Cambridge University Press, 2003.
22. A. L. Mackay, Acta Crystallogr., 1962, 15, 916.
23. L. D. Marks, Philos. Mag. A, 1984, 49, 81.
24. J. E. Jones, Proc. R. Soc. London, Ser. A, 1924, 106, 463.
25. R. P. Gupta, Phys. Rev. B: Condens. Matter, 1981, 23, 6265.
26. F. Cleri and V. Rosato, Phys. Rev. B: Condens. Matter, 1993, 48, 22; J. P. K. Doye, D. J. Wales and M. A. Miller, J. Chem. Phys., 1998, 109, 8143; J. P. K. Doye and D. J. Wales, Phys. Rev. Lett., 1998, 80, 1357.
27. J. P. K. Doye and F. Calvo, Phys. Rev. Lett., 2001, 86, 3570; J. P. K. Doye and F. Calvo, J. Chem. Phys., 2002, 116, 8307.
28. E. Panizon and R. Ferrando, Phys. Rev. B: Condens. Matter, 2015, 92, 205417.
29. V. A. Sharapov, D. Meluzzi and V. A. Mandelshtam, Phys. Rev. Lett., 2007, 98, 105701.
30. W. Schottky, Phys. Z., 1922, 23, 448; E. F. Westrum, J. Chem. Thermodyn., 1983, 15, 305.
31. P. H. E. Meijer, J. H. Colwell and B. P. Shah, Am. J. Physiol., 1973, 41, 332.
32. L. Rubinovich, M. I. Haftel, N. Bernstein and M. Polak, Phys. Rev. B: Condens. Matter, 2006, 74, 035405.
33. F. Sciortino, W. Kob and P. Tartaglia, J. Phys.: Condens. Matter, 2000, 12, 6525; T. V. Bogdan, D. J. Wales and F. Calvo, J. Chem. Phys., 2006, 124, 044102; D. J. Wales, Chem. Phys. Lett., 2013, 584, 1.
34. J. J. Moré, B. S. Garbow and K. E. Hillstrom, User guide for MINPACK-1, Technical Report ANL-80-74, Argonne Nat. Lab., Argonne, IL, 1980; E. Jones, T. Oliphant and P. Peterson, et al., SciPy: Open source scientific tools for Python, 2001.
35. V. Vitek and T. Egami, Phys. Status Solidi B, 1987, 144, 145.
36. B. Zhu and M. Hou, Eur. Phys. J. D, 2012, 66, 1; K. Laasonen, E. Panizon, D. Bochicchio and R. Ferrando, J. Phys. Chem. C, 2013, 117, 26405.
37. D. Schebarchov and D. J. Wales, Phys. Chem. Chem. Phys., 2015, 17, 28331.
38. Z. Li and H. A. Scheraga, Proc. Natl. Acad. Sci. U. S. A., 1987, 84, 6611.
39. D. J. Wales and J. P. K. Doye, J. Phys. Chem. A, 1997, 101, 5111.
40. B. E. Husic, Dopant effects on solid-solid transitions in atomic clusters, M. Phil., University of Cambridge, 2015.
41. V. A. Mandelshtam and P. A. Frantsuzov, J. Chem. Phys., 2006, 124, 204511.
42. V. A. Sharapov and V. A. Mandelshtam, J. Phys. Chem. A, 2007, 111, 10284.
43. V. A. Mandelshtam, P. A. Frantsuzov and F. Calvo, J. Phys. Chem. A, 2006, 110, 5326.
44. D. Schebarchov and S. C. Hendy, Phys. Rev. Lett., 2005, 95, 116101; D. Schebarchov and S. C. Hendy, Phys. Rev. B: Condens. Matter, 2006, 73, 121402; D. Schebarchov and S. C. Hendy, Eur. Phys. J. D, 2007, 43, 11.
45. D. Schebarchov and D. J. Wales, Phys. Rev. Lett., 2014, 113, 156102.
46. M. Bixon and J. Jortner, J. Chem. Phys., 1989, 91, 1631.
47. A. B. Harris, J. Phys. Chem., 1974, 7, 1671; Y. Imry and M. Wortis, Phys. Rev. B: Condens. Matter, 1979, 19, 3580.
48. M. E. Fisher and A. E. Ferdinand, Phys. Rev. Lett., 1967, 19, 169; Y. Imry, Phys. Rev. B: Condens. Matter, 1980, 21, 2042.
49. J. D. Honeycutt and H. C. Andersen, J. Phys. Chem., 1987, 91, 4950.
50. A. Sachdev, R. I. Masel and J. B. Adams, Z. Phys. D: At., Mol. Clusters, 1993, 26, 310; J. P. K. Doye, Phys. Rev. B: Condens. Matter, 2003, 68, 195418.
51. I. L. Garzón and A. Posada-Amarillas, Phys. Rev. B: Condens. Matter, 1996, 54, 11796.
52. 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.
53. K. Michaelian, N. Rendón and I. L. Garzón, Phys. Rev. B: Condens. Matter, 1999, 60, 2000.
54. C. Kittel, Introduction to solid state physics, Wiley, 2005.
55. D. Faken and H. Jónsson, Comput. Mater. Sci., 1994, 2, 279.
56. V. Rosato, M. Guillope and B. Legrand, Philos. Mag. A, 1989, 59, 321.
57. E. Aprà, F. Baletto, R. Ferrando and A. Fortunelli, Phys. Rev. Lett., 2004, 93, 065502.
58. F. Calvo, J. P. K. Doye and D. J. Wales, Nanoscale, 2012, 4, 1085.
59. K. Sutherland-Cash, D. J. Wales and D. Chakrabarti, Chem. Phys. Lett., 2015, 625, 1; F. Calvo, D. Schebarchov and D. J. Wales, J. Chem. Theory Comput., 2016, 12, 902.
60. J. P. K. Doye and D. J. Wales, J. Chem. Phys., 1995, 102, 9659.

---

Footnote

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

---