Grand and Semigrand Canonical Basin-Hopping.

We introduce grand and semigrand canonical global optimization approaches using basin-hopping with an acceptance criterion based on the local contribution of each potential energy minimum to the (semi)grand potential. The method is tested using local harmonic vibrational densities of states for atomic clusters as a function of temperature and chemical potential. The predicted global minima switch from dissociated states to clusters for larger values of the chemical potential and lower temperatures, in agreement with the predictions of a model fitted to heat capacity data for selected clusters. Semigrand canonical optimization allows us to identify particularly stable compositions in multicomponent nanoalloys as a function of increasing temperature, whereas the grand canonical potential can produce a useful survey of favorable structures as a byproduct of the global optimization search.

Introduction

Structure prediction is essential in many areas of computational science, ranging from molecular physics and biochemistry to soft and condensed matter. For a given system with definite size, the global optimization problem is usually nontrivial owing to high-dimensional potential energy landscapes, and many methods have been proposed to locate low-energy configurations.[1−4] The global minimum is fundamentally important and often carries essential insight into the interactions responsible for the emergence of specific morphologies, and it plays an important role in explaining self-assembling motifs and symmetries.[5,6] However, in many applications, temperature can play a significant role and the entropic contribution to the free energy of individual configurations becomes important. Examples of entropy-driven structural transitions have been reported in atomic clusters,[7,8] proteins,[9] colloids,[10] glasses,[11] and pressurized materials.[12] The determination of configurations that are low in free energy can proceed by the a posteriori analysis of molecular simulations, often employing biases in order to sample the energy landscape more efficiently and based on system-dependent order parameters.[13−17] Such free energies are global and can encompass many potential energy minima, as expected in the context of phase transitions. Local free energies can also be defined for individual isomers, which are related to global quantities through suitable grouping procedures that require additional knowledge about the connectivity of the minima.[18,19]

It is possible to locate the global free energy minimum among an existing database of structures by evaluating the entropy using the harmonic approximation[20,21] and, when affordable, incorporating anharmonic corrections.[21−23] More recently, calculating the free energy directly on-the-fly during global optimization was proposed,[24] producing a promising procedure for exploring energy landscapes, since the free energy minimum was encountered faster this way than by postprocessing a sample based on optimizing the potential energy alone. In the present contribution, we further extend this approach by addressing systems with variable size or composition, which should be treated in the grand canonical or semigrand canonical ensembles, respectively. Fluctuations in the number of particles occur in the case of nucleation of fluids and their absorption into porous materials, as well as in the increasingly important problem of reversible gas storage for energy production. These ensembles correspond to situations in which the system exchanges particles with a (possibly fictitious) reservoir, thereby controlling size or composition at fixed temperature.

Grand canonical ensembles are characterized by the chemical potential (or chemical potential difference) and a finite temperature, with the Gibbs free energy being the potential of interest that, in turn, controls the size or composition around equilibrium. As the chemical potential varies, changes in the Gibbs free energy are indicative of different regimes in which the system grows or shrinks or reaches equilibrium values in absorption isotherms. In semigrand canonical ensembles, particularly stable compositions should be manifested by plateaux in the segregation isotherms. Compared to free energy global optimization, the need to calculate local Gibbs free energies for systems with varying size or composition requires sampling these additional variables as well. The extra degrees of freedom further justify the use of the harmonic approximation to approximate the entropy component in a computationally efficient manner. In practice, the harmonic approximation requires calculating the vibrational frequencies at local minima, which involves constructing and diagonalizing the dynamical matrix (the mass-weighted Hessian).

Our results for atomic clusters in the grand canonical ensemble indicate that increasingly large clusters are obtained as the chemical potential is increased or the temperature is decreased, in agreement with conventional nucleation theories. The results of grand canonical basin-hopping simulations in the harmonic approximation are also found to agree with a model for the grand canonical partition function fitted to reproduce the size-dependent heat capacity for specific clusters. Our other example application deals with model nanoalloys treated in the semigrand canonical ensemble, which is at fixed total size but varying alloy composition. Here, we show that the semigrand canonical basin-hopping method efficiently locates the stability plateaux in the composition isotherms that were previously reported based on alternative simulation methods.[25]

The article is organized as follows. The next section describes the method in its general formulation and details the harmonic expression employed for the local Gibbs free energy associated with individual potential energy minima. Practical details regarding the implementation of the basin-hopping method are also given in relation to Monte Carlo moves that change the system size. The results for the grand canonical ensemble are presented in Section and discussed in the context of the homogeneous nucleation problem. The semigrand canonical application to model nanoalloys is described in Section , followed by concluding remarks in Section .

Methods

Grand Canonical Formulation

In the grand canonical ensemble, the volume , temperature T, and chemical potential μ are fixed, whereas the pressure, energy, and number of particles, N, can fluctuate. The grand partition function describing an equilibrium distribution is thenwhere β = 1/kBT, with kB the Boltzmann constant, and is the canonical partition function of the system with fixed number of particles, N. In previous work, we have constructed as a sum over contributions from local minima, i.e., the superposition approach.[20,21,26−29] Here, the classical vibrational density of states can be written aswhere κ(N) = 3N – 6 is the number of nonzero eigenvalues for the Hessian matrix, v̅α = [∏κ(να(j)]1/κ( is the geometric mean vibrational frequency of minimum α, να(j) is the normal-mode frequency of the jth mode in this minimum, and Eα is the corresponding potential energy. The approximation in eq corresponds to using harmonic vibrational frequencies. The superposition approach can also incorporate quantum effects[23] and anharmonicity,[21−23,28,30−32] but it is usually employed in the harmonic approximation to obtain a rapid survey of thermodynamic properties, which is guaranteed to be ergodic by construction.

In recent work, we have demonstrated how the superposition framework can usefully be applied within grand and semigrand canonical formulations to examine equilibrium thermodynamics.[25] In the present contribution, we show how a potential function based on the grand (and semigrand) canonical ensembles can be used in the context of basin-hopping global optimization. This approach is a natural extension of the free energy basin-hopping method,[24] generalized so that the size (or composition) is permitted to change. To find the largest contribution to the grand potential from local minima of any size N, we adapt the acceptance criterion to use the potentialwhere we have included the rotational partition function for a rigid rotor with inertia tensor I, for completeness. A translational contribution was not included since we are considering nontranslating clusters. The sign definition in eq enables us to formulate the location of a maximum contribution to the grand potential as a global minimization. In a monatomic system, nα( = 2N!/oα, where oα is the order of the point group.[26,33−35] In this harmonic/rigid rotor approximation, the occupation probability p for structures containing N atoms can then be written in terms of the sum over the corresponding subset of minima (denoted by )

A Metropolis acceptance criterion was applied at each step, using an acceptance probability based on min{1,exp[−(ξα(new) – ξα(old))/kBTGCBH]}, where TGCBH is a fictitious temperature parameter that determines how often uphill moves in ξ are accepted. Along with the temperature parameter, another key choice in basin-hopping that affects efficiency is the coordinate perturbation scheme applied before each local minimization. Here, we employed perhaps the simplest scheme based on Cartesian coordinate displacements, drawn from a uniform distribution with a fixed maximum value. Many other possibilities have been considered in the literature, along with variations in the acceptance condition, and the present approach could be combined with any of these methods in future work. So long as the key local minimization is included,[2,36,37] efficiency gains might be possible. Here, we adopted one other modification, since moves that involve changes in N are likely to be much more disruptive than geometrical perturbations. We therefore considered moves changing the number of atoms only at intervals of Δ basin-hopping steps. Before changing the number of atoms, the structure of the current minimum in a Markov chain over blocks was saved, along with the value of ξα. After adding or removing an atom, the structure obtained after minimization was used as the initial seed for a local Markov chain of Δ coordinate perturbations, accepting moves according to the condition min{1,exp[ – (ξα(new) – ξα(old))/kBTGCBH]} for fixed N. At the last step of each block, the minimum with the lowest value of ξα was used in a block accept/reject test and compared with the current minimum in the block Markov chain. If the block move was rejected, then the current structure was reset to the one saved in the block Markov chain before proposing another move that changed the number of atoms. The same BH temperature parameter was employed for both types of acceptance check, although different values, or indeed different criteria, could certainly be considered.

Two different schemes were compared for grand canonical basin-hopping (GCBH) in the present work to analyze the role of different contributions to the potential. The first scheme used ξα, as defined above, and the second scheme based the sampling on ξα(0) = Eα – Nμ, omitting the last term of the right-hand side of eq involving the canonical partition function of the minimum. In effect, this version corresponds to neglecting entropic contributions. The results are compared in Section .

The grand canonical potential in eq can be reinterpreted as a semigrand canonical potential for a binary system with variable numbers of A- and B-type particles, NA and NB, respectively, but N = NA + NB fixed. The potential employed in the accept/reject criterion is thenwhere Δμ = μB – μA is the chemical potential difference, and nα( = 2NA!NB!/oα. In Section , we describe a semigrand canonical basin-hopping (SGCBH) scheme with and without the rotational contribution (), geometric perturbations disabled, and particle insertion/deletion moves replaced by exchange moves that transmute particles.

Steps That Change the Cluster Size

Steps within the blocks of constant size were proposed in a manner similar to that in previous work, including procedures to move weakly bound or surface atoms.[37] Steps to change size were proposed by adding or deleting single atoms, with probability p+ and p– = 1 – p+, respectively. All of the results reported below simply used p+ = p– = 0.5 throughout. Atomic clusters bound by a pair potential were chosen for the first application of grand canonical basin-hopping. For pair potentials (but also for embedded-atom potentials, as used in our semigrand canonical basin-hopping examples), we can easily identify the most weakly bound atom for each minimum, and this was the atom removed in steps that reduced the cluster size. To add an atom, the center of coordinates was first located, along with the largest atomic radial distance, rmax. An atom was then added at a random point on the sphere with radius rmax + δr (where δr = σ for Lennard-Jones systems) by generating a three-component vector with each entry drawn from the normal distribution with zero mean and unit variance[38] and then normalizing appropriately. Gaussian random variables were generated using the Box–Muller algorithm.[39]

Following each atom addition or deletion, the resulting configuration was immediately minimized; if this quench failed, then the attempted size change was simply rejected. Local minima were also rejected if they did not correspond to connected single clusters. Here, we used a depth first search[40] to check for a percolating network of atoms, which proved to be particularly useful in previous studies of clusters bound by short-range potentials.[41] The present results therefore exclude fragmented systems, focusing on single clusters.

Application of Basin-Hopping to a Grand Canonical Potential

We considered clusters bound by the Lennard-Jones (LJ) potential,[42] where the potential energy isϵ is the equilibrium pair well depth, 21/6σ is the equilibrium pair separation, and r is the distance between particles i and j. Reduced units with ϵ = σ = 1 are used throughout.

Depending on the values of T and μ, the GCBH runs either exhibit sizes that shrink to a single atom or grow quite steadily to contain hundreds of atoms for the run lengths considered here. To prevent indefinite growth, an upper limit of 1000 atoms was applied.

Model for the Size-Dependent Density of States

In order to interpret the observed behavior, a model for the size-dependent partition function was constructed for comparison with the GCBH runs using data available for selected LJ clusters. Here, the input included results for predicted global minimum energies, point groups, normal modes, and moments of inertia up to LJ1610.[37,43−46] We also employed the heat capacity curves calculated by parallel tempering Monte Carlo[47,48] for N = 13, 31, 55, 75, 129, 135, and 309. The objective was to predict relative values for the canonical partition function as a function of N, and we represented the N-dependent terms of interest aswhere the subscript 0 refers to the putative global minimum for N atoms, and Stirling’s approximation was used for ln N!. The first term is the usual Boltzmann factor, associated with the global potential energy minimum for N. It serves as the lower bound to the potential energy distribution of local minima at each size. The next two terms derive from the harmonic approximation to the vibrational density of states, again for the global minimum, which serves as a baseline. The fourth term is the N-dependent part of the rotational partition function for the global minimum, and the N! term comes from the number of distinct permutation-inversion isomers, with the corresponding point group order in the final term. This last term is a representation of the quadrature over the potential energy density of local minima and their associated vibrational and rotational densities of states. The quadrature has H slices, where a geometric progression with common ratio R(N, T) = exp{[γ(N) – Δ(N)/kBT]/H} has been summed. It subsumes variation of the rotational and vibrational densities of states relative to the reference values for the global minimum, which are divided out in eq . The contribution of the global minimum is included separately as 1/o0, since this structure is likely to have particularly high symmetry,[6,26] and the spectrum of energy minima is better represented as discrete rather than continuous at low energy. Δ(N) is the potential energy range spanned by minima containing N atoms, and the parameter γ(N) includes the rate of increase in the number of minima with potential energy. Hence, ΔE(N)/H is the width of each quadrature slice, and the contribution of slice q in the canonical partition function is exp{[γ(N) – ΔE(N)/kBT]q/H}.

In selecting the above model, we have been guided by the densities of minima calculated in previous work[31] for LJ31 and LJ75. The energy range and growth rate, ΔE(N) and γ(N), were fitted to best reproduce the heat capacity curves from parallel tempering Monte Carlo data, using ΔE(N) = ΔE0 + NΔE1 + N2/3ΔE2 and γ(N) = γ0 + Nγ1 + N2/3γ2. Various other representations were also considered; the choice in eq produces an acceptable fit for the melting peak in the heat capacity curves (Table ). The fitting is only qualitative, since we make no attempt to reproduce additional peaks due to phase-like transformations below the melting point.[26,29,37,49−57] We require only a qualitative description of the growth in the number of local minima, convoluted with the change in rotational and vibrational densities of states relative to the global minimum reference values. The results presented below correspond to H = 30 quadrature slices, which is sufficient to converge the properties of interest. The resulting heat capacities for some selected cluster sizes are shown in Figure .

Table 1

Fitted Parameters Employed for the Model Canonical Partition Functions

ΔE0	8.34177	γ0	28.7118
ΔE1	0.332144	γ1	0.56546
ΔE2	–0.452931	γ2	0.47651

Figure 1

Variation of the heat capacity, C/N, with temperature for some selected cluster sizes using the model defined by eq . The peak shifts to higher temperature with increasing size, and the results illustrated correspond to N = 13, 31, 38, 55, 75, 129, 309, and 1610.

It is instructive to consider the individual summation terms in eq associated with a particular N at selected fixed values of the chemical potential μ and temperature T. In Figure , we plotwith evaluated from the model defined in eq and representing a relative (i.e., not normalized) probability of cluster size N for a given μ, , and T. Note that test runs where the rotational contribution was omitted gave essentially the same results (data not shown).

Figure 2

Cluster size distributions, as quantified by defined in eq , for (a) fixed μ = −7.5 and temperatures of 0.05, 0.075, 0.1, 0.15, 0.2, 0.25, and 0.3 (reduced units) and (b) fixed T = 0.05 and chemical potentials −7.1, – 7.2, – 7.3, – 7.4, – 7.5, and −7.6.

The nonmonotonic trends in Figure exhibit a minimum, which corresponds to the least probable cluster size Nlp, i.e., the N that minimizes . (The relationship between Nlp and the definition of a critical nucleus will be investigated in more detail in future work.) Figure shows that, for a range of μ and T values, Nlp(μ, T) is consistently intermediate between Nmin = 3 and Nmax = 1610. In contrast, Figure shows that the value of N that maximizes , denoted by Nmp as short-hand for “most probable”, always appears either at the maximum (Nmax) or the minimum (Nmin) size considered. Since Nmp accounts for the largest contribution to and always occurs at one end of the size range considered, one result of a GCBH run is a prediction of cluster growth or dissociation for a given μ and T. If a survey of low-lying minima is not required, then each run can be terminated once the growth or dissociation is clear.

Figure 3

Least probable cluster size, Nlp, as a function of μ and T. For the largest and smallest values of μ considered, Nlp corresponds (or is very close) to one limit of the cluster size range. However, for intermediate values of μ, it can occur elsewhere, corresponding to the minima in Figure .

Figure 4

Most probable cluster size, Nmp, as a function of μ and T. Note that Nmp always occurs at (or very close to) one limit of the cluster size range considered.

GCBH Results

We find that cluster growth is favored by low temperature and larger μ, and the values of μ above which GCBH runs produced growth, μ*(T), are summarized as a function of temperature in Table . The same values were obtained in runs of 107 GCBH steps with size moves attempted every 50 or 1000 steps. The corresponding block sizes for accept/reject testing after a size move were set to 40 and 900. Runs that exhibited cluster growth did not reach the limit of 1000 atoms. The maximum Cartesian coordinate perturbation was fixed at 0.4, and the temperature parameter employed in both the standard BH steps (constant size) and in the block accept/reject tests was fixed at TGCBH = 1.5 reduced units throughout in this initial survey. The temperature parameter must be carefully distinguished from the physical temperature that appears in the grand potential. As in standard basin-hopping, effective values for the temperature parameter should be large enough for the algorithm to escape from low-lying local minima and small enough for low-lying regions to be properly explored. The GMIN code includes a replica exchange option where exchanges are attempted between parallel runs with different values of the basin-hopping temperature parameter,[58] but we have not explored this extension in the present work.

Table 2

For Fixed Temperatures, the Value of μ Above Which GCBH Runs Produced Cluster Growth, μ*(T)a

T	μ*(T)
0.001	–5
0.1	–3
0.15	–2
0.2	–2
0.4	0

The μ values considered were in integer steps.

The observation of a minimum in the relative probability size distribution for ranges of T and μ seems to be consistent with the appearance of a critical nucleus in classical nucleation theory.[59−62] In the latter framework, clusters containing fewer particles than the critical value are predicted to shrink and clusters containing more particles should exhibit spontaneous growth. However, we note that various subtleties exist in the treatment of cluster nucleation by simulation,[63,64] and connections with the procedure adopted in the present work require further investigation. Our results simply represent a search to optimize the potential ξα, defined in eq , with no attempt to describe a vapor/droplet equilibrium. For reference, recent calculations of the chemical potential for saturation (equilibrium), μ, with alternative potentials for argon produce values around −8.7 and −10.4 for kBT/ϵ ≈ 0.61 and 0.78, respectively.[65]

Runs using ξα(0) = Eα – Nμ, which do not require normal-mode analysis, gave results consistent with the low-temperature limit of ξα, with μ* = −5. The lowest minima located as a function of cluster size are an interesting and potentially useful byproduct of the GCBH runs. Since systematic global optimization is not being performed explicitly for each cluster size, we would not expect to locate the true global minima for all N in the range encountered. Nevertheless, searches that involve changes in size, and the block move structure, may provide a useful survey of the low-lying morphologies. To illustrate this possibility, the lowest potential energies encountered in two selected GCBH runs compared to the lowest known minima are plotted as a function of N in Figure . These plots include a range of μ values for temperatures of 0.001 and 0.1 and illustrate the trends that also appear for higher temperatures. At T = 0.001, runs with μ values large enough for cluster formation often encounter the lowest known minimum up to N around 200 atoms. For larger sizes, the potential energy difference generally increases with N.

Figure 5

Potential energy difference (ΔV) in ϵ between the lowest minimum encountered in selected GCBH runs and the lowest known minimum as a function of cluster size N. Top: T = 0.001; bottom: T = 0.1. The common μ values for the two panels are −2 (red), −3 (orange), −4 (yellow), −5 (green), and −6 (blue). The black line in the lower panel corresponds to μ = −1. In both cases, the results are for GCBH runs of 106 steps with size moves attempted every 1000 steps, block size for accept/reject testing of 900 steps, maximum Cartesian coordinate perturbation fixed at 0.4, and TGCBH = 1.5 throughout (reduced units). The data points are joined to guide the eye.

For T = 0.1, the lowest minima encountered are close to the global minimum for smaller sizes, and the energy difference again generally rises with N. However, for μ = −2, we see the difference decrease above about N = 370. These trends basically reflect the proportion of steps that the GCBH run spends around each cluster size. The better a given N value is sampled, the smaller the deviation from the energy of the putative global minimum. The results for higher temperatures (omitted for brevity) are consistent with this pattern. To obtain a survey of low-energy structures as a function of size, we would therefore choose a low temperature, together with a large enough value of μ to produce clustering.

Application of Basin-Hopping to a Semigrand Canonical Potential

Model and Calculation Details

To demonstrate semigrand canonical basin-hopping (SGCBH) and compare with previous work,[25] we consider icosahedral AgPd55– nanoalloys modeled by the same many-body (Gupta) potential (with the same parameter values). Using independent SGCBH runs, each at a different value of Δμ, we seek the composition and mixing pattern that minimize the value of ξα defined in eq , with the rotational contribution omitted, as in ref (25), or included for comparison. Before proceeding, we note here that the calculations reported in ref (25) employed a different definition of ξ based on Δμ(NAg – NPd) = Δμ(NAg – (N – NAg)) = 2ΔμNAg – NΔμ. The chemical potential difference used in this reference should, therefore, be multiplied by two for comparison with the present results.

Each SGCBH run consisted of 105 basin-hopping steps at kBTSGCBH = 0.01 eV, with the stoichiometry reset to NAg = 28 and species distribution randomized every 500 steps to enhance exploration.

Equilibrium Composition

The Ag content of the lowest encountered configuration obtained at the three physical temperatures of 0, 100, and 300 K is plotted versus Δμ in Figure . At low temperature, the segregation isotherms display particularly striking plateaux for NAg = 12 and 42, as well as secondary steps at intervening sizes. These plateaux indicate particularly stable compositions and are associated here with highly symmetric motifs in which the nanoalloy is enriched in silver at the surface, either at the vertices (NAg = 12) or entirely (NAg = 42). Those motifs were also identified in ref (25), and after scaling the chemical potential difference by a factor of 2, the results agree well with those obtained here.

Figure 6

Ag content in the lowest lying AgPd55– icosahedron at a given (relative) chemical potential Δμ = μAg – μPd, found using SGCBH based on the potential defined in eq for three temperatures. In the ball-and-stick representation of the most persistent stoichiometries, Ag atoms are gray (lighter) and Pd atoms are magenta (darker).

In ref (25), the harmonic superposition method was applied to the same system, and the equilibrium fraction of silver atoms contains all of the important thermal fluctuations. Among the contributions to the grand partition function included in the superposition method, the semigrand canonical basin-hopping method seeks only the term corresponding to the potential energy minimum that makes the largest contribution. Only at T = 0, when thermal occupation of higher energy minima is completely suppressed, do the two calculations coincide precisely. We show in Figure the importance of such fluctuations by comparing the SGCBH and HSA results at 300 K, correcting the previous results by scaling the chemical potential difference, as explained above. In this graph, we have also superimposed the basin-hopping results obtained when the rotational contribution to the partition function is included.

Figure 7

Ag content in the lowest lying AgPd55– icosahedron at a given (relative) chemical potential Δμ = μAg – μPd, found using semigrand canonical basin-hopping (SGCBH) based on the potential defined in eq for T = 300 K, and compared to the results of harmonic superposition approximation (HSA) that include the contribution of a database of minima (solid line). The results of SGCBH accounting for the rotational partition function are also shown.

The neglect of higher energy minima in the SGCBH method explains the staircase character of the silver fraction in Figures and 7, whereas the superposition method gives smoother variations even at 100 K. As temperature increases, thermal fluctuations also increase and further smoothen the equilibrium silver fraction.[25] A similar effect is found for the SGCBH approach, which results from the increasingly similar Gibbs free energies of the various isomers. While the fluctuations predicted by the HSA appear relatively significant at 300 K, they are actually well represented by the global minimum of the semigrand partition function. It is also worth noting that some of the differences between the SCGBH and HSA calculations are related to the inclusion of nonicosahedral structures in the latter results, which contribute significantly to the stabilization of the composition 43:12 in Ag/Pd near Δμ = 1.2 eV.[25] This discrepancy observed between the superposition and basin-hopping calculations would likely be reduced if the SGCBH simulations were no longer restricted to sampling the icosahedral homotop.

Finally, Figure confirms that the rotational partition functions plays essentially no role on the segregation isotherms for this cluster. This result was expected because all minima are homotops of the two-layer Mackay icosahedron; hence, the moments of inertia differ only in the ways the atomic masses are allocated within this framework.

Some Remarks on the High Symmetry of Stable Compositions

The high symmetry of the most persistent compositions is consistent with the principle of maximum symmetry,[6,26,66] which suggests that structures with a higher symmetry content are most likely to lie in the low- and high-energy tails of the distribution. Since exceptions to this principle are not difficult to find, we do not expect it to be universally applicable, but there seems to be no doubt that many global minima observed for a wide variety of systems have high symmetry content. An explanation can be found by writing the total energy as a sum over contributions from a many-body expansion, involving single atom, pairwise and three-body terms, etc. If these terms are drawn from the same distribution, then geometrical symmetry (degeneracies) would be manifested as correlation. The variance is larger when correlation is present. Symmetrical structures are therefore more likely to have particularly high or particularly low energy. Low-lying structures are therefore likely to exhibit symmetry.

More formally, denote the mean and variance of a variable, X, drawn from probability distribution p(X) as μ and σ2. Contributions to the total energy are correlated when symmetry is present, and for the corresponding probability distributions, p(X,X′) ≠ p(X) p(X′). If we sum M terms with mean μ and variance σ2 from the same distribution, then the variance of the mean is Mσ2 + M(M − 1)ρσ2 where the correlation ρ is defined byFor ρ = 0, the variance is Mσ2, but for ρ = 1, it rises to M2σ2. On average, we therefore expect to find structures with a higher symmetry content in the tails of the distribution.[6]

We have previously exploited symmetry-biased moves in the context of basin-hopping global optimization, and the core orbits scheme produced efficiency gains of 1 or 2 orders of magnitude for some benchmark atomic clusters.[67] The present results suggest that the maximum symmetry principle may also apply to nanoalloys.

Conclusions

Grand canonical ensembles are commonly encountered in problems relevant to nucleation or absorption. In this article, a methodology is introduced to search directly for the configuration that minimizes the Gibbs free energy and therefore makes the largest individual contribution to the grand potential. The method is based on basin-hopping global optimization and extends a recent effort[24] where the free energy landscape for a given system is minimized using the harmonic approximation for the entropy at finite temperature. Here, we have further developed the method to treat a variable number of particles (grand canonical ensemble) or, at fixed size, the composition in heterogeneous systems (semigrand canonical ensemble). Two applications illustrating both situations were presented, dealing with the nucleation of LJ clusters at fixed temperature or with the progressive segregation in nanoalloys as the chemical potential is varied.

In the grand canonical case, the method predicts that the global minimum in the local Gibbs free energy switches from dissociated states to bound clusters for larger chemical potentials and lower temperatures. Those results are in agreement with a specific model for the size- and temperature-dependent canonical partition function fitted to reproduce the heat capacities of reference clusters. They are also consistent with one of the primary conclusions of nucleation theory, namely, the existence of a critical nucleus that minimizes the Gibbs free energy.[59−62]

Our semigrand canonical application to silver–palladium nanoalloys can be compared with the predictions of the harmonic superposition method,[25] which is similar to the present approach but incorporates the contributions of different minima to the grand partition function, not just the largest individual terms. The two approaches generally agree well, with discrepancies appearing at higher temperature, where new families of minima are stabilized, which are omitted in the present work.

Compared with the superposition method for the same statistical ensembles, grand canonical basin-hopping is less demanding because only the global minimum of the corresponding ensemble is sought instead of entire (ergodic) samples covering all relevant regions of the free energy landscape for the different system sizes or compositions. For the semigrand canonical ensemble, our approach provides a powerful way of solving the combinatorial problem of chemical ordering for a nonrigid lattice at finite temperature.

Future applications of the present methodology include the important case of absorption into porous materials, especially at low temperature and high densities, for which conventional grand canonical simulations are not practical. In particular, the contribution of zero-point motion can be straightforwardly incorporated in the expression of the partition functions from knowledge of the individual vibrational frequencies. Moreover, the notorious difficulty of sampling high-density phases would be alleviated owing to the systematic local minimization step of the basin-hopping procedure, thus increasing the chance of accepting the corresponding moves.