Metallic nanoparticles meet MetadynamicsL. Pavan, K. Rossi and F. Baletto1, a)

Physics Department, King’s College London, WC2R 2LS, UK

(Dated: July 2015; Revised October 2015)

We show how standard Metadynamics coupled with classical Molecular Dynamics can be successfully appliedto sample the configuration space of metallic and bimetallic nanoclusters. This is achieved via the implemen-tation of collective variables related to the pair distance distribution function of the nanoparticle. As paradig-matic examples we show an application of our methodology to Ag147, Pt147 and their alloy AgshellPtcore at2:1 and 1:1 chemical compositions. The proposed scheme is able to reproduce the known solid-solid structuraltransformation pathways, based on the Lipscomb’s diamond-square-diamond mechanisms, both in the monoand bimetallic nanoparticles. An estimate of the free energy barrier for those mechanisms is provided anddiscussed.

I. INTRODUCTION

Mono and bimetallic nanoparticles (NPs) find a widenumber of applications ranging from catalysis andbiomedicine to optoelectronics and magnetic data stor-age due to their high surface to volume ratio, pecu-liar shape and d-band shift.1,2 Nanoclusters’ chemophys-ical properties strongly depend on the interplay betweentheir size, morphology, and chemical composition. Un-derstanding the thermal stability of a configuration, in-cluding complex entropic contribution difficult to addressexperimentally,3 is highly desiderable and many theoret-ical attempts have been presented in the literature.4–6

The investigation of the solid-solid structural pathwaysand the evaluation of the free energy barriers (FEBs)among different configurations may shed light on howcluster chemical features vary due to ageing and externalfactors such as temperature or pressure. Although therecent developments made possible to detect solid-solidtransformations via electron beam irradiation7,8, insightson the atomistic details and an accurate sampling of theenergy landscape can be provided only using a numericalapproach9.

Numerical technique tools to obtain a quantitativesampling of the nanoclusters’ conformational space canbe divided in two families. Double ended searches arebased upon the foreknowledge of the initial and finalpoint of the transition and consist in producing the min-imum energy path, using an eigenvector-following algo-rithm to get the transition state and then a steepest-descent energy minimisation to achieve the lowest energymechanism. Notable examples are the nudged elasticband,10 and the double-ended transition path sampling(DETPS).11 Double ended approaches have been appliedto a wide range of systems which also include bimetallicnanoclusters12. These calculations, however, become ex-pensive when used to analyse structural rearrangementsin systems with high friction, high free energy barriers,and/or a very rough free energy landscape13. Temper-ature accelerated and biased sampling techniques repre-

a)Electronic mail: francesca.baletto@kcl.ac.uk

sent an opposite approach: an initial configuration of thesystem is excited or perturbed and forced to visit new andunknown isomers. Metadynamics14,15, adaptive biasingforce16, umbrella sampling17 and parallel tempering18 areall renown techniques based upon this idea. Perturba-tion methods have been commonly used to detect order-disorder transition while rarely for simulate solid-solidstructural transformations in nanoalloys.19–24

In this paper, we will show how Metadynamics(MetaD) can be successfully employed to sample the con-figuration space of relatively large metallic and bimetal-lic nanoparticles at room temperature. MetaD algo-rithm coarses the system dynamics in a collective vari-ables (CVs) space where a history dependent poten-tial is exploited to accellerate rare events as minimum-minimum transitions and to reconstruct the system’s freeenergy surface (FES) projection in this order parameterspace.14,15 It has been shown that this method can be ap-plied to finite inorganic nanosystems such as semiconduc-tor/quantum dots25,26, alkali halides nanostructures27,28

Lennard-Jones29,30. However it has been applied only tometallic nanoclusters with less than 20 atoms31,32 andnever to bimetallic cases. Our MetaD scheme, with CVsbeing window functions on the pair distance distribu-tion function reproduces the well known five and sixdiamond-square-diamond (DSD)33,34 mechanism in bothmono and bimetallic systems. Moreover Mackay’s pre-diction of the transformation of an icosahedron into ananticuboctahedron has been proven numerically for thefirst time.

The section ”Model and Method” reports the method-ology adopted to explore the NP free energy landscape;”Result and Discussion” contains first a focus on the de-scritions of the structural pathways detected during ourMetaD runs, followed by a quantitative analysis of theirenergy barriers. A resume of our results and discussionis presented in the final section.

II. MODELS AND METHOD

We consider monometallic (Ag,Pt) and bimetallic(AgPt) nanoparticles at 2:1 and 1:1 chemical composi-

2

FIG. 1. Closed shell polyhedra at 147 atoms. From leftto right, top row: icosahedron, decahedron, cuboctahedron.Bottom row: hexagonal closed packed and anticuboctahe-dron.

tion, made of 147 atoms. At this size, noble and quasi-noble metals are likely to assume closed shell polyhe-dra geometries. The more common are the icosahedron(Ih), decahedron (Dh), cuboctahedron (Co), while theless are the anticuboctahedron (aCo) and common hexag-onal close-packed geometries (hcp). These NP architec-tures are all reported in Figure 1. Ag and Pt present asmall size mismatch (less than 4%) and they show a largeimmiscibility gap in the bulk. Taking into account thatAg has a lower surface energy than Pt, an AgshellPtcoreordering is expected35.

We perform classical molecular dynamics (MD) sim-ulations, coupled with MetaD, where a velocity Verletalgorithm is used to solve Newton’s equations of motion.A time step of 5 fs is considered, the temperature is con-trolled via an Andersen thermostat. Atomic interactionsare modelled within the second moment approximationof the tight binding theory36 (TBSMA), where the po-tential VTBSMA(x) is the sum of atomic contribution,Ei

TBSMA:

EiTBSMA =

nv∑j 6=i

Aabe−pab

(rij

r0ab

−1)

(1)

−

√√√√ nv∑j 6=i

ξ2abe−2qab

(rij

r0ab

−1)

.

Here the sum is extended up to the number of atomsnv within an appropriate cut-off distance from atom i;a and b refer to the chemical species of the two atoms,r0ab is the bulk nearest neighbour distance for the ho-mometallic atomic pair and their arithmetical averagefor the bimetallic case. Aab and ξab are related to thecohesive energy, while pab and qab and their product de-termine the range of the repulsive and of the attractivepart of the potential. Details of their parametrizationand their value for Ag and Pt are reported in Ref.37 whilefor, AgPt parameters are taken from Ref.35.

MetaD enables the enhanced sampling of the free en-ergy surface of a system. It couples a coarse graineddescription of the system in the CV space to an history-dependent biasing potential, ∆V , which evolves as thesum of Gaussians of height ω and width σ added everytG time interval

V (x) = VTBSMA(x) + ∆V (S(x), t) (2)

∆V (S(x), t) =∑

t′=tG,2tG..

ωe−[S(x)−s(t′)]2

2σ2 .

S(x) represents the set of collective variables chosen anddefines the order parameter space where the MetaD po-tential evolves; s(t′) is the value of the collective variableat time t′. The efficiency and physical faithfulness ofMetaD is strongly related to the choice of a sensible setof CV. The metric described by the chosen CVs should beable to distinguish between the various configurations thesystem visits during a run. The CV should also be rep-resentative of all slow degree of motions involved in thevarious structural transformations of interest. Further-more their number should be limited in order to explorea low dimensionality space thus avoiding undesired com-putational expenses. If one of the above conditions is notrespected the system may be forced to visit high energyregions and/or conformations and processes of interestmay be hidden.15 When the MetaD potential reaches adiffusive regime, the FES projection in the CV space canbe reconstructed as the negative of the MetaD potential∆V (S(x), t) in Eq3.

The research of a good order parameter able to dif-ferentiate among the various isomers of a cluster is nottrivial38,39. We have noted that an almost complete in-formation on NP morphology is encoded in the pair dis-tance distribution function (PDDF) of the nanoclusters,as depicted in Figure 2. We choose our CV set accordingto the following criterion: it must be able to distinguishthe geometries in Fig.? as these are the main structuresof interest in our research. Thus, we introduce specifi-cally tailored order parameters corresponding to windowfunction (WF ) on the PDDF constructed via a sigmoidfunction:

WF (x) =∑

i,j;i 6=j

1−(

rij−d0

r0

)61−

(rij−d0

r0

)12 , (3)

where rij is the distance between atom i and atom j,r0 the window width, always set to 0.05 times the bulklattice constant and d0 is the window function charac-teristic distance. d0 is set to be 1.354 and 3.4 of thebulk lattice parameter respectively for the stacking faultnumber (SFN) and the maximum pair distance difference(MPDD), as highlighted in Figure 2.

A qualitative physical argument can be provided tosupport our choice. The SFN is related to a character-istic hcp peak in the PDDF. This is a topological defect

3

FIG. 2. Pair distance distribution function for icosahedral(orange), decahedral (blue), cuboctahedral (grey), and hcp(black) motifs of a 147 atoms cluster. The selected windowsaround the stacking fault number (SFN), and the maximumpair distribution difference (MPDD), are highlighted in yellowand grey shadowing, respectively.

obtained due to the intersection of two planes with dif-ferent symmetry orientation. At the nanoscale, as in thebulk, phase transformations happen via stacking faults,as demonstrated for the diamond-square-diamond mech-anism in small Lennard Jones clusters40. Consequentlya bias on this CV should not constrain the FES explo-ration to unphysical structural transition or highly ener-getic states and should represent the sliding and rotationof (111) planes: the slow degree of freedom in the DSDmechanisms. The MPDD is set where the pair distribu-tion function of the different geometries is less overlap-ping. In such a way we are confident of being able todiscriminate at least the main structural morphologies.

We resort to the use of these two window functions onspecific lattice distances because they adopt different val-ues for each of the geometries of interest in our researchand are computationally very cheap. This will enableus to easily extend our methodology to larger systems? .We note that a small loss in accuracy in our calcula-tion may be caused by the partial correlation of the two.We remark that this set of CVs can be applied with apartial confidence to bimetallic nanoalloys: the differentchemical species are not treated explicitly, hence struc-tural transformations involving chemical reordering arerarely reproduced biasing on those two reaction coordi-nates. As we focus on a stable AgshellPtcore ordering,SFN and MPDD will capture the structural transitionpathways of interest.

NP geometrical features are monitored on-the-fly viacommon neighbour analysis (CNA)41 during the courseof our MetaD + MD simulations. The CNA associates tothe local network of each nearest neighbours atom paira signature made of three integers (r, s, t). r is the num-ber of common nearest neighbour, s the number of bondsbetween the r-common nearest neighbours and t is thelongest chain among the s-bonds. Different CNA signa-

FIG. 3. Common neighbour network (yellow) of the red-coloured atomic pair considered. Left panel shows a (422)signature, associated with twin boundaries, while the rightpanel depicts a (555) signature, related to pairs lying on a5-fold symmetry axis.

tures are able to distinguish whether a pair of atoms isin a bulk environment (s ≥ 4) or on the surface (s=2,3);whether the bulk is crystallographic(e.g. a FCC bulk hasonly (4,2,1)); if the pair belongs to a fivefold axis ((5,5,5)signature). Any isomer is characterised by a typical per-centage value of a few signatures. The polyhedra consid-ered in Figure 2 can be distinguished in terms of their(4,2,2) and the (5,5,5) signatures, illustrated in Figure 3.At 147 atoms, the (4,2,2) % decreases from 39% (Ih) to25% (Dh and hcp) and then zero (Co); while the (5,5,5)% ranges between 5.2 % (Ih), 0.9 % (Dh) and zero (Coand hcp).

A significant change in the CV and in the CNA signa-tures is therefore distinctive of structural transition, asin the paradigmatic example of an Agshell92Ptcore55 runsketched in Figure 4, where three different basins, Ih (or-ange), Co (blue) and Dh (pink), are explored. Structuresinvolved in any presumed structural rearrangements arethen quenched in order to identify precisely the initial,saddle, and final configurations. The lower panel of Fig-ure 4 reports the quantity ∆Equench defined as the po-tential energy of a structure, rescaled with respect to thetotal energy of the Ih minimum, which results to be themost energetically favorable in all systems. Potential en-ergy barriers (PEB) are then calculated by means of theDEPTS available in the OPTIM package,11,42 where therelaxed structures obtained during a MetaD run, serveas the initial and final configurations.

III. RESULTS AND DISCUSSION

All the results shown in the following have been ob-tained at 300K. Gaussians 0.25 eV high are depositedwith a period of 20ps. Their width σ is 15 along the SFNdimension and 10 for the MPDD. Such values are basedon the standard deviations of the CV during unbiasedmolecular dynamics simulations at 300 K and are chosento insure an efficient conformational flooding.

4

FIG. 4. Agshell92Ptcore55 MD + MetaD simulation outcomeand analysis. From top to bottom: Total energy difference(∆Etot = Etot − Equench

Ih ), CV values, on-the-fly CNA signa-tures, Total energy difference of relaxed structures taken fromthe simulation (∆Equench = Equench − Equench

Ih ). A colourscheme is applied to discriminate the independent conforma-tions basin explored: Co-basin in blue, Ih-basin in orange,Dh-basin in pink.

A. Solid-Solid Structural Transformation Pathways

The calculated solid-solid transformation mechanismsfrom Ih to Dh, Co and aCo basins ,and viceversa, arereported in Figures 5, 6 and 7, respectively. The mo-tion of atoms in a (111) facets is clearly identified via thedifferent coloring while the initial, saddle and final con-figurations are shown. As expected, they correspond to aDSD mechanism as predicted by Mackay33. Further, wehighlight that the MetaD transition pathways obtainedfor monometallic systems is identical to the one obtainedvia OPTIM/DETPS.

A DSD mechanism consists of a stretching and rota-tion of triangular facets into a diamond first and then asquare by a collective screw dislocation of atoms. Thedislocation corresponds to a rotation of different anglesaccording to the initial and the final configuration.

FIG. 5. Initial, saddle and final configurations of DSD mech-anism for Ih into Dh transformation. Multicolored atoms de-limit a facet in the original Ih. Rotation along the depictedaxis is of 36 degrees.

FIG. 6. Initial, saddle and final configurations of DSD mech-anism for Ih into Co transformation. Multicolored atoms de-limit a facet in the original Ih. Rotation along the depictedaxis is of 60 degrees.

Conversely a squared facet can transform into two tri-angular facets by the opposite movement. Five parallelo-grams are involved in the collective rearrangement Dh↔Ih (Figure 5) and six in the case of the Co ↔ Ih (Figure6) and of an aCo ↔ Ih transformation (Figure 7).

An aCo can convert into an Ih throughout a rotationin opposite ways of the two triangular facets perpendic-ular to the aCo three-fold rotation axis, by 60 degreeswith reference to each other, along the same axis. Thedifference with the six DSD mechanism between Co↔ Ihis that now two opposite parallel triangular facets rotateabout their normal and three pairs of abutting triangularfacets remain unchanged, rotating about the axis that isperpendicular to their common edge and belongs to thetwin plane. According to our knowledge, this is the first

FIG. 7. Ih transforming into aCo via DSD mechanism. Saddlepoint is pictured in between the two. Multicolored atomsdelimit a facet in the original icosahedron which turns intoa diamond one by rotating by 60 in opposite direction withrespect to the parallel triangular facet.

5

FIG. 8. Initial, saddle and final configuration in the case ofDh into Ih transformation in Ag92shellPt55core. Ag atoms arein grey, Pt in blue. Atoms from Ag shell are removed to showPt inner core.

FIG. 9. Initial, saddle and final configuration in the case ofDh into Ih transformation in Ag74shellPt73. Ag atoms arein grey, Pt in blue. Red lines highlight facets evolution instructural rearrangement process.

time Mackay’s path has been found in a numerical simu-lation.

The solid-solid transition pathways for perfect core-shell ordering are the same as described earlier. More-over, it is important to note that the transition happenssimultaneously in the two chemical species, as shown inFigure 8. No other mechanism appears to be able toconnect the Ih and Co, Ih and aCo and Ih to Dh. ForAgshell74Pt73, a few transitions as a Dh into Ih, and anIh into a defected Co have been reproduced. We notethat in this system structural rearrangements involvesboth DSD and surface reconstruction mechanisms as de-picted Figure 9. Further studies are then needed in orderto fully understand the interplay between the collectivegeometrical mechanism and atomic inter-diffusion.

In Ag147, the Ih, Co and the aCo basins explorationreaches a diffusive regime, we also observed the Dh con-version into an Ih, but the opposite transition has notbeen visible. In the case of the Pt147 clusters only thetransformations into a Ih geometry from Dh Co, andaCo basins have been observed before the explorationof highly defected structures. For Agshell92Ptcore55, thetransformation Co into Ih, and Dh into Ih have been de-tected together with their inverse processes. On the otherhand, only the conversion from the aCo into the Ih basinhas been seen.

B. Energy Landscapes Analysis

In this subsection we will discuss the reconstruction ofthe free energy landscapes for the various systems. Signif-icance of entropic effects at T=300K will also be adressedvia a comparison between free energy and potential en-ergy barriers, when possible. We dedicate a paragraphto each system considered.

Ag147

For the Co → Ih (Co ← Ih) interconversion the FEBis of 0.4eV (3.2eV). In the case of aCo → Ih (aCo → Ih)transformation, it is of 0.5eV (3.0eV). The value of thesame processes evaluated via DEPTS result respectivelyof 0.5eV (3.2eV) for the Co→Ih and of 0.64eV (3.49eV)for the aCo→ Ih structural rearrangements.

In the case of the Dh to Ih conversion, the inverse struc-tural transition has not been observed. We introduced anon-canonical approach to obtain an (over)estimate ofthe DhrightarrowIh FEB. In this case, the non con-vergence is likely caused by the fact that defected Ihstructures, e.g. rosette-like43 overlap and hide the Dhbasin. However, the quantities in our on the fly geomet-rical analysis, CNA and ∆Equench, can be used as ad-ditional order parameters to identify precisely when theFEB around the Dh basin have been overcome and theIh basin exploration has started. The free energy land-scape may be then reconstructed as the inverse of theMetaD potential needed to exit from the initial Dh basinto explore the Ih one. We approximate the free energybarrier to the highest energy point in that landscape.Following this scheme, the Dh → Ih free energy barrieris (over)estimated to be of 0.6 ev. Error on our overes-timate may be heuristically taken to be at least equal tothe height of one Gaussian (0.25 eV). Activation barrierfor the Dh →Ih transition calculated via DEPTS are of0.45 eV and 2.17 eV for the direct and inverse processrespectively. In the case of the fully converged Ih ↔ Coand Ih ↔ aCo the entropic contribution calculated asthe difference between FEB and PEB is of the order of0.1 eV. . If we consider a similar entropic effect on thisprocess the calculated upper bound and relative error forthe estimate of that free energy barrier results sensible.

Pt147

For this system only rearrangements from Dh, Co andaCo to Ih have been simulated. From Ih, only defectedgeometries have been detected, this is because the pres-ence of highly distorted icosahedral motif is even morelikely in Pt than in Ag. The simulation results also sug-gest that in Pt147 FEB from Ih to defected structures arelower with respect to the one for perfect geometries as Coand Dh. The system remains trapped into the large Ihconformational basin and convergence cannot be reached

6

in a reasonable time. Our CVs in fact do not hinderthe transformation in defected geometries which mightoverlap the Dh or Co basin. We recall that Pt147 showsa strong tendency to form Ih defected structure, relatedto its ability to reconstruct facets as described in Ref43,which appear to be entropically favourite even with re-spect to the perfect Ih structure at high temperatures.This may be the physical origin that prevents the systemto reach a diffusive regime between the perfect geometryconformations.

FEBs have been estimated within the non-conventionalapproach described in theAg147 section. They result of1.9 eV for Dh → Ih, 1.7 eV for Co→Ih and 2.07 eV foraCo→Ih, while the corresponding PEB values are of 2.11eV, 1.9 eV and 1.87 eV respectively. Notwithstanding thenon conventional approach FEB barriers for Co-Ih andDh-Ih transition result slightly lower than PEB ones, in asimilar fashion to the one calculated for the converged Agruns. If we take in account that the structural rearrange-ment pathways from MetaD and DETPS algorithms arealmost identical, and further assume that the entropiceffects in Ag and Pt systems may be comparable, thismakes us confident that the CVs employed are actuallysensible and that our free energy barrier estimate is rea-sonable.

AgPt147

In the case of the Co↔Ih transformation inAgshell92Ptcore55 the calculated FEB is 0.25eV and 3.6eV respectively. The details of this MetaD run are re-ported in Figure 4 and the FES reconstruction is shown inFigure 10. Roughness of the landscape can be bestowedto the sampling of intermediate defected structures. Asmall barrier for escaping the Co basin may be under-stand as this shape is energetically very unfavorable.

The Dh into Ih runs are also convergent but a partialbasins overlap appear due to the exploration of a widenumber of defected structures thus causing a significantunderestimation of the free energy barrier if the canonicalapproach is used. Hence we again adopt the FEB upperbound evaluation previously described. The upper boundof the Dh → Ih (Dh ← Ih) transition is of 1.5 eV (4.3eV).

FEB investigation for Ag74shellPt73 is discouraged bythe non diffusivity of the simulation, by the fact that thestructural rearrangement pathway is more complex withrespect to the other system and by the simple fact thatthey will depend on the arbitrarily chosen chemical ar-rangement of the initial structures used in the simulation.

The energy barriers for the Dh-Ih core-shell structuralrearrangement processes is in between the ones found forpure Ag and pure Pt clusters. Co-Ih and aCo-Ih barriersin Agshell92Ptcore55 appear smaller than the one calcu-lated for monometallic systems as these configurationsare essentially not energetically stable.

FIG. 10. Co and Ih basins free energy surface reconstructionfor an Ag92shellPt55core. Hot colours in the CV space high-light an high free energy of the corresponding configuration.

IV. CONCLUSION

We have shown how Metadynamics can be adapted tothe study of solid-solid structural transitions in metallicand bimetallic NPs by introducing specific CV. These aretwo window functions set onto the pair distance distri-bution function, where they refer to the stacking faultnumber and the maximum pair distribution distance.

We have shown that the solid-solid structural tran-sitions among the most common closed shell polyhe-dra happen following the Lipscomb’s diamond-square-diamond mechanism, which is a collective screw dislo-cation motion. Further, we have demonstrated that theMackay’s description of the interconversion of an Ih intoan aCo throughout a DSD mechanism is reproducible.According to our in-silico experiment structural rear-rangement mechanisms that allow to connect minima onthe free energy surface are identical to the one on the po-tential energy landscape. We would like to remark thatthe DSD mechanism appears to be a universal pathwayas it takes place in monometallic as well as in nanoalloyswith a small either non negligible mismatch, such as inAgAu44 and AgPt, this work.

Entropic effects appear to play a minor role at roomtemperature as the calculated free energy barriers arecomparable in magnitude to the activation energy ones.We acknowledge that the overestimate of our methodol-ogy may be compensating for some of the entropic effectswhich would reduce further the magnitude of the FEB.We notice that Ih and Dh basin may overlap in the con-sidered CV space due to the exploration of defected Ihstructures which may hidden the Dh basin. Finally, weremark that our set of collective variables may be ap-plied to monometallic (e.g. Ag, Pt) and bimetallic sys-tems with a small mismatch (e.g. AgPt) and/or witha well defined core-shell chemical ordered. However, in

7

nanoalloys with a large mismatch a limitation occurs: abroadening of the peaks in the PDDF may result in theimpossibility of uniquely defining the characteristic dis-tances in Eq. 3.

ACKNOWLEDGMENTS

LP and FB, KR thank the financial support byU.K. research council EPSRC, under Grants No.EP/GO03146/1 and EP/J010812/1. KR acknowledgesfinancial support by U.K. research council EPSRC,Grants No. ER/M506357/1. The simulations were car-ried out using the Faculty and Departmental computa-tional facilities at Kings College London.

1F. R., J. J., and R. L. Johnston, Chem. Rev. 108, 845 (2008).2F. Baletto and R. Ferrando, Rev. Mod. Phys. 77, 371 (2005).3D. J. Wales and P. Salamon, Proc. Natl. Acad. Sci. U. S. A.111, 617 (2014).

4A. S. Barnard and P. Zapol, J. Chem. Phys. 121, 4276 (2004).5R. S. Berry and B. M. Smirnov, Physics-Uspekhi 48, 345(2005).

6Z. H. Li and D. G. Truhlar, Chem. Sci. 5, 2605 (2014).7Z. W. Wang and R. E. Palmer, Phys. Rev. Lett. 108, 245502(2012).

8C.-C. Chen, C. Zhu, E. R. White, C.-Y. Chiu, M. C. Scott,B. C. Regan, L. D. Marks, Y. Huang, and J. Miao, Nature496, 74 (2013).

9E. Ringe, R. P. V. Duyne, and L. D. Marks, J. Phys. Chem.C 117, 15859 (2013).

10G. Henkelman and H. Jo, J. Chem. Phys. 113, 9901 (2000).11D. J. Wales, Molecular Physics 102, 891 (2004).12F. Calvo, A. Fortunelli, F. Negreiros, and D. J. Wales, J.Chem. Phys. 139, 111102 (2013).

13P. G. Bolhuis, D. Chandler, C. Dellago, and P. L. Geissler,Ann. Rev. Phys. Chem. 53, 291 (2002).

14A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U. S. A.99, 12562 (2002).

15A. Laio and F. L. Gervasio, Rep. Prog. Phys. 71, 126601 (2008).16E. Darve, D. Rodrıguez-Gomez, and A. Pohorille, J. Chem.Phys. 128, 144120 (2008).

17G. Torrie and J. Valleau, J. Comput. Phys. 23, 187 (1977).

18D. J. Earl and M. W. Deem, Phys. Chem. Chem. Phys. 7, 3910(2005).

19G. Adjanor, M. Athenes, and F. Calvo, Europ. Phys. J. B53, 47 (2006).

20J. P. Neirotti, F. Calvo, D. L. Freeman, and J. D. Doll, J.Chem. Phys. 112, 10340 (2000).

21F. Calvo and C. Mottet, Phys. Rev. B 84, 035409 (2011).22F. Calvo, J. P. K. Doye, and D. J. Wales, Nanoscale 4, 1085(2012).

23F. Calvo and E. Yurtsever, J. Chem. Phys. 140, 214301(2014).

24E. G. Noya and J. P. K. Doye, J. Chem. Phys. 124, 104503(2006).

25C. Bealing, R. Martonak, and C. Molteni, J. Chem. Phys.130, 124712 (2009).

26C. Bealing, R. Martoak, and C. Molteni, Solid State Sciences12, 157 (2010).

27D. Bochicchio and et al., J. Chem. Phys. 140, 2024911(2014).

28J. Liu and Q. Ge, J. Phys. Cond. Mat. 23, 345401 (2011).29F. Pietrucci and W. Andreoni, Phys. Rev. Lett. 107, 85504(2011).

30G. A. Tribello, J. Cuny, H. Eshet, and M. Parrinello, J.Chem. Phys. 135, 114109 (2011).

31G. Santarossa, A. Vargas, M. Iannuzzi, and A. Baiker, Phys.Rev. B 81, 174205 (2010).

32L. Pavan, C. Di Paola, and F. Baletto, Eur. Phys. J. D 67,24 (2013).

33A. L. Mackay, Acta Cryst. 15, 916 (1962).34W. N. Lipscomb, Science 153, 373 (1966).35L. O. Paz-Borbon, R. L. Johnston, G. Barcaro, and A. For-

tunelli, J. Chem. Phys. 128, 134517 (2008).36V. Rosato, M. Guillope, and B. Legrand, Philos. Mag. A 59,321 (1989).

37F. Cleri and V. Rosato, Phys. Rev. B 48, 22 (1993).38D. J. Wales, J. Chem. Phys. 142, 130901 (2015).39F. Pietrucci and R. Martoak, J. Chem. Phys. 142, 104704(2015).

40J. Uppenbrink and D. J. Wales, Faraday Trans. 87, 215 (1991).41J. D. Honeycutt and H. C. Andersen, J. Phys. Chem. 91,4950 (1987).

42S. Trygubenko and D. Wales, J. Chem. Phys. 120, 2082(2004).

43E. Apra, F. Baletto, R. Ferrando, and A. Fortunelli, Phys.Rev. Lett. 93, 65502 (2004).

44C. F. and R. L. Johnston, Appl. Phys. Lett. 92, 023112 (2008).