Size and shape-dependent melting mechanism of Pd nanoparticles
Movaffaq Kateb, Maryam Azadeh, Pirooz Marashi, Snorri Ingvarsson
NNoname manuscript No. (will be inserted by the editor)
Size and shape-dependent melting mechanism of Pdnanoparticles
Movaffaq Kateb · Maryam Azadeh · Pirooz Marashi · Snorri Ingvarsson
Received: date / Accepted: date
Abstract
Molecular dynamics simulation is employed to understand the ther-modynamic behavior of cuboctahedron (cub) and icosahedron (ico) nanopar-ticles with 2 −
20 number of full shells. The original embedded atom method(EAM) was compared to the more recent highly optimized version as inter-atomic potential. The thermal stability of clusters were probed using potentialenergy and specific heat capacity as well as structure analysis by radial distri-bution function ( G ( r )) and common neighbor analysis (CNA), simultaneously,to make a comprehensive picture of the solid state and melting transitions. Theresult shows ico is the only stable shape of small clusters (Pd – Pd usingoriginal EAM and Pd using optimized version) those are melting uniformlydue to their small diameter. An exception is cub Pd modeled via optimizedEAM that transforms to ico at elevated temperatures. A similar cub to icotransition was predicted by original EAM for Pd – Pd clusters whilefor the larger clusters both cub and ico are stable up to the melting point.As detected by G ( r ) and CNA, moderate and large cub clusters were showingsurface melting by nucleation of the liquid phase at (100) planes and growth ofliquid phase at the surface before inward growth. While diagonal (one cornerto another) melting was dominating over ico clusters owing to their partitioned M. KatebScience Institute, University of Iceland, Dunhaga 3, IS-107 Reykjavik, IcelandE-mail: [email protected]. AzadehDepartment of Mining and Metallurgical Engineering, Amirkabir University of Technology,Tehran, Iran
Present address: School of Metallurgy and Materials Engineering, Universityof Tehran, Iran
P. MarashiDepartment of Mining and Metallurgical Engineering, Amirkabir University of Technology,Tehran, IranS. IngvarssonScience Institute, University of Iceland, Dunhaga 3, IS-107 Reykjavik, Iceland a r X i v : . [ c ond - m a t . m t r l - s c i ] F e b Movaffaq Kateb et al. structure which retarded the growth of the liquid phase. The large ico cluster,using optimized EAM, presented a combination of surface and diagonal melt-ing due to the simultaneous diagonal melting started from different corners.Finally, the melting temperature as well as latent heat of fusion were calcu-lated and compared with available models and previous studies which showed,unlike the present result, the models failed to predict size-dependent motifcrossover.
Keywords
Size-dependent · Nanoparticle · Melting · Enthalpy
Palladium nanoparticles have received growing interest due to their novel prop-erties such as size (Chen et al, 2016) and shape (facet) (Wang et al, 2017) de-pendent catalytic activity, controlled oligomerization (Zhivonitko et al, 2016)and enhanced hydrogen storage (Rangel et al, 2016). However, due to thehigher surface to volume ratio (Schmidt et al, 1998; Safaei et al, 2008), thenanoparticles are known as the most unstable structures among different nano-solids. Thus, it is necessary to predict their thermal stability for practicalapplications specifically when instability at elevated temperature can be con-sidered as a failure. For instance, reshaping, melting or aggregation of catalystclusters might happen far below the bulk melting point ( T mb ) which results inchanging their properties as well as functionality.The nanoparticles melting temperature ( T mp ) and enthalpy ( H mp ), amongthe other thermodynamic properties, have received considerable attention (Qi,2016). Several models have been developed to predict the nanoparticles size-dependent T mp (Goldstein et al, 1992; Jiang et al, 1999; Safaei, 2010) and H mp (Zhang et al, 2000; Jiang et al, 2002; Attarian Shandiz and Safaei, 2008;Fu et al, 2017). However, developing a universal model requires an extensiveeffort and each model has to be verified with experimental or simulation data(Liang et al, 2017) which, to our knowledge, is not available for Pd clusters. Inaddition, these models fail to describe the melting dynamic and size-dependentmelting mechanisms. For instance, quite recently, it has been shown that re-shaping of Ag particles depends on the size and a transition from homogeneousto surface melting mechanism occurs by an increase in the particle size (Lianget al, 2017).Alternatively, atomistic Monte Carlo (MC) (Westergren and Nordholm,2003) and molecular dynamics (MD) (Baletto et al, 2002; Pan et al, 2005; Miaoet al, 2005; Schebarchov and Hendy, 2006) simulations have proven to be anexcellent tool for probing the nanoparticles stability and melting behavior. Forinstance, Baletto et al (2002) have shown that for small Pd cluster icosahedron(ico) structure is more stable while decahedron (dec) and cuboctahedron (cub)are more stable at the moderate and large clusters, respectively. In this regard,they compared Rosato-Guillop´e-Legrand and embedded atom method (EAM)force fields those were in close agreement. However, they only focused on theenergetics of the different structures rather than showing a dynamic transition ize and shape-dependent melting mechanism of Pd nanoparticles 3 between allotropes. Thus, it still remains unclear that under what conditionan unstable shape transforms into a more stable counterpart. Later, Pan et al(2005) reported in the case of Pd cluster both of the cub and ico shapesare stable in a wide range of temperatures using Sutton-Chen (SC) potential.They also reported surface melting mechanism which was detected using radialdistribution function, G ( r ). However, the G ( r ) has been defined separately foreach shell of the cluster in their work which can be misleading i.e. the diffusionof the surface atoms, those have more degree of freedom, into inner shells causesa broadening of the G ( r ) peak for the outer shell which can be misinterpretedas surface melting (Bertoldi et al, 2017). They regretted this drawback byreporting quasi-solid core at the melting point ( T mp ) of clusters while their G ( r ) patterns indicated broad peaks for all shells (Pan et al, 2005). They alsoobserved a minor peak prior to the main melting peak in the specific heatcapacity C p of the ico cluster which they attributed to the surface melting.Using a similar potential, Miao et al (2005) reported the surface melting for aspherical Pd cluster based on bond-orientational order parameters (BOP).The BOP is a local structure characterization method obtained by determiningthe angle between a list of nearest neighbors (NN) (Steinhardt et al, 1983).Thus, it determines which atoms are in solid/liquid states and enables thestudy of the surface melting more precisely. Later, Schebarchov and Hendy(2006) pointed out the BOP is unable to distinguish between fcc and hcp(Steinhardt et al, 1983) and it is not efficient to detect solid state transitionsin the clusters. To maintain brevity, cub to ico transition can be detectedby mapping twining boundaries in ico cluster those can be characterized ashcp atoms. Thus, they suggested common neighbor analysis (CNA), insteadof BOP, which is very sensitive to the angles between different pairs of NN(Tsuzuki et al, 2007). The efficiency of CNA for the structural transition hasbeen demonstrated earlier by Baletto et al (2000, 2001) and interested reader isreferred to Ref. (Rossi and Ferrando, 2007) to learn about the potential of CNAfor the detailed characterization of clusters. They successfully demonstrateddec to ico transition through a solid-liquid state for a Pd cluster. This isinteresting since they used the EAM force field and obtained a different resultthan static energy calculation by Baletto et al (2002) with similar potential. Itis worth noting that for small clusters those have a high degree of freedom andmight undergo several structure transformations at elevated temperatures, theLindemann index (Lindemann, 1910) ( δ L ) gives scattered results as reportedelsewhere (Alavi and Thompson, 2006; Zhang and Douglas, 2013).In the present work, the phenomenological melting of symmetric shell Pdclusters of cub and ico structure is investigated. Since EAM potential predictssolid-state phase change, according to previous studies, original and highly op-timized EAM were utilized in the present study. In addition, different methodsof caloric curves, G ( r ) and CNA were utilized to make a comprehensive un-derstanding of the nanoparticles melting. This accompanied with a brief dis-cussion on the limitation of each method and the resulting misinterpretation.Finally, the size-dependent results were compared to state of the art models. Movaffaq Kateb et al. (Plimpton, 1995; Plimpton andThompson, 2012).The original (Foiles et al, 1986) and highly optimized (Sheng et al, 2011)EAM was used to depict the inter-atomic potential between Pd atoms. Eq. (1)represents the original formulation of EAM potential: E i = F i (cid:88) i (cid:54) = j ρ ij ( r ij ) + 12 (cid:88) i (cid:54) = j U ij ( r ij ) (1)where E i and F i are cohesive and embedding energies of atom i , respectively. ρ ij ( r ij ) is the electron density of j atoms located around the i atom at dis-tance r ij . Clearly, F i is a many-body interaction term while U ij takes the pairinteraction into account.The EAM potential was confirmed for describing solid characteristics suchas cohesive energy and elastic constant as well as metals melting point (Foileset al, 1986; Sheng et al, 2011). Moreover, it is reliable in determining thetransitional properties especially the heat of fusion and heat capacities aboveroom-temperature. The EAM has also been verified for a quantitatively cor-rect description of such nanoscale systems; for instance, surface energy andgeometry of low index surfaces.The time integration of the equation of the motion was performed regardingthe Verlet algorithm (Verlet, 1967; Kateb and Dehghani, 2012) with a timestepof 3 fs. The temperature control was done using the Nose-Hoover thermostatwith a damping of 30 fs. These conditions designed to generate positions andvelocities sampled from canonical (NVT) ensemble. The initial velocities of theatoms were defined randomly from a Gaussian distribution at the appropriatetemperature of 300 K and relaxed for 300 ps in NVT ensemble. Since practicalnanoparticle melting is mostly performed in the vacuum, the heat associatedwith the particle’s melting cannot be removed so efficiently by the surroundingmedium. But the particles are in contact with the substrate which allowscontrol over temperature and hence, the NVT ensemble provides a realisticrepresentation of such systems. The simulations were performed by startingat 300 K, and then the temperature was elevated at a heating rate of 1 . × K/s.2.2 Cluster preparationThe Pd nanoparticles were considered to be in the cub and ico forms thosefound in experimental characterizations (Jos´e-Yacam´an et al, 2001). For dif- version 14 Apr 2013 available at http://lammps.sandia.govize and shape-dependent melting mechanism of Pd nanoparticles 5 ferent sizes of given clusters, cub and ico were made considering the so-calledmagic number (Poole Jr and Owens, 2003) or a total number of cluster atoms( N t ) which is described as the function of the shell number ( n ): N t = 13 (10 n + 15 n + 11 n + 3) (2)where n = 0 denotes a mono-atomic system and n ≥ ∼ n = 2 −
20 shells( N t = 55 − n = 8)and 8-ico as an example. The figure also shows a slice of particles indicating8-cub is made by fcc and surface atoms while the 8-ico presents more com-plicated arrangement i.e. 20 fcc tetrahedral (green) those surrounded an icoatom (indicated by yellow in Fig. 1d) separated by twinning grain boundariesdetected as hcp (red).We chose to use the cluster radius ( R p ) defined by Guinier formula (Miaoet al, 2005) which is nearly constant during surface diffusion, solid-state tran-sition and at the beginning of melting. R p = R g (cid:114)
35 + r a (3)where R g is particle gyration radius and r a is atomic radius equal to 0.137 nm.The term particle size hereafter is referred to as D = 2 R p calculated at ambienttemperature.2.3 Melting criteriaThe melting is commonly detected from isotherm transition in the so-calledcaloric curve(s) e.g. potential energy ( U ) plotted versus temperature ( T ) (Shimet al, 2002; Zhao et al, 2001). Qi et al (2001) defined T mp as the temperaturewith the maximum apparent heat capacity. In MD simulations, specific heatcapacity at constant pressure ( C p ) can be calculated by derivative of the U average ( U ave ) during heating: C p = dU ave dT + 32 (cid:60) (4)here (cid:60) stands for the universal gas constant.It is worth noting that the precision of T mp value calculated from C p isdependent on the number available data points ( U vs. T ) and a limited numberof sampling, which is usually the case in MC simulation, results in a huge error.Several methods are used in MD simulations to identify the melting pro-cess based on atomic specifics. The first criterion proposed by Lindemann(Lindemann, 1910) who stated the melting of crystals occurs when the aver-age amplitude of atomic vibrations is higher than a threshold value. The global δ L is a system average of atomic quantity which shows a linear increase withtemperature increment in solid-state and a step change due to the melting. Movaffaq Kateb et al.
A BC D
Fig. 1
Illustration of (a) 8-cub and (b) 8-ico particles with 8 being the shell number, n .The crystal arrangement is shown through a slice of (c) 8-cub and (d) 8-ico indicating fcc,hcp, ico and unknown using green, red, yellow and blue colors, respectively. However, most of the vibrations of the surface atoms in the small clusters,which have more degree of freedom, were assumed as melting behavior by thismodel (Alavi and Thompson, 2006; Zhang and Douglas, 2013). This is a seri-ous issue since it may lead to misinterpreting of the surface melting instead ofa solid state transition. The most straightforward structure analysis is offeredby G ( r ) or pair correlation function (Iida and Guthrie, 1988). It describes howdensity varies as a function of distance in a system of particles from a referenceparticle. G ( r ) = (cid:104) πr ρ a dr (cid:105) T (5)where ρ a is the atom numbers density, r is the distance from reference particleand dr determines the bin size. The angle brackets i.e. (cid:104)(cid:105) T denote the timeaverage at constant T . ize and shape-dependent melting mechanism of Pd nanoparticles 7 Here G ( r ) was extracted by 300 ps relaxation of clusters in NVT ensembleat desired T i.e. 300 – 2500 K with 100 K steps. Then the G ( r ) averagedout over the whole time with 10 − nm bin size without periodic boundarycondition (PBC).This results in a pattern of several peaks corresponding to number anddistance of NNs which applies to a wide range of materials. The melting tran-sition causes a variation in the density and can be detected by shifting andbroadening of peaks in the G ( r ) pattern. However complex solid-state transi-tion at nanoscale such as cub to ico with constant coordination number andeven distance, is very hard to determine with G ( r ).Recently, CNA has shown promising tool due to providing the possibilityof a distinction between allotropic transitions and melting process. Briefly,the CNA determining local crystal structure based on the decomposition of1st NNs obtained from G ( r ) into different angles (Faken and J´onsson, 1994;Tsuzuki et al, 2007). It is quite sensitive technique to the symmetry of differentpairs of bonds. Thus a twining grain boundary as the main difference of icoand cub cluster can be determined based on a slight angle difference between1st NNs while it holds entire properties of an fcc atom.2.4 VisualizationThe open visualization tool (OVITO) package was used to generate atomisticillustrations (Stukowski, 2009). ∆E , is optimized to be equal to 0.02 eV/atomin agreement with thermodynamically assessed value at room temperature(Dinsdale, 1991) and more precise than 0.026 eV/atom obtained by originalEAM (Foiles et al, 1986). This value is very important for realization of icoto cub transition and stability of ico structure. It is believed the vacancy for-mation energy, E v , is the elementary mechanism of melting and thus it is oneof the most important parameter to capture realized melting. The optimizedforce field results in vacancy formation energy of 1.4 eV similar to value esti-mated from experimental melting point (Kraftmakher and Strelkov, 1970) andbetter than 1.44 eV obtained from original potential (Foiles et al, 1986). Thevalues for surface energy, γ sv , in the optimized EAM predicts higher valuesthan original EAM regarding all plains. The experimental value of the γ sv isestimated from contact angle representing an average over all planes, but it is Version 2.7.1 available at http://ovito.org/ Movaffaq Kateb et al.
Table 1
The values obtained by optimized EAM (Sheng et al, 2011) force field in com-parison with experimental or ab initio results and original EAM(Foiles et al, 1986). The a denotes lattice parameter.Parameter experiment EAM EAM(unit) or ab initio original optimized T mb (K) 1828 (Rao and Rao, 1964) 1680 1750 a ( nm ) 0.389 (Ashcroft and Mermin, 1976) 0.389 0.389(100) 2000 (Tyson and Miller, 1977) 1370 1645 γ sv (mJ/m ) (110) 2000 1490 1747(111) 2000 1220 1529 ∆E (eV/atom) fcc-hcp 0.02 (Dinsdale, 1991) 0.026 0.02 E v (eV/atom) 1.4 1.44 1.4 still higher than both EAM values. Finally, the optimized EAM results in T mb value closer to experimental value (Rao and Rao, 1964).Fig. 2 demonstrates the result of both EAM potential in determiningthe cohesive energy, E c , of solid Pd in comparison with ab initio results.The minimum in fcc denotes equilibrium a and E c equal to 0.389 nm and3.911 eV/atom, respectively, which both EAM are in agreement with ab ini-tio result. However, the a and E c values for hcp determined by ab initio,are smaller than both EAM by 0.005 nm and 0.06 eV, respectively. In general,compared to ab initio results, both EAM are slightly off in the repulsive range,but they are acceptable in attraction part. The optimized EAM, however, tendto zero more smoothly.Further, the dynamic structure factor, S ( q ), was calculated which is aconventional method to examine liquid structure (Iida and Guthrie, 1988).The S ( q ) was calculated here using the Fourier transform of G ( r ). S ( q ) = 1 + ρ a (cid:90) ( G ( r ) − e iqr dr (6)here q is the wave-vector in reciprocal space.The G ( r ) was obtained in the same way described in subsection 2.3 butwith dr = 10 − nm and PBC for 500 atoms in total. The wave-vector incrementwas chosen to be compatible with PBC i.e. 2 π/L with L being the cube sidelength equal to 5 × a (1.945 nm).Fig. 3 shows S ( q ) for the liquid Pd calculated using both EAM in compari-son with tight binding calculation (Alemany et al, 1999) and the experimentalresult obtained by x-ray Raman scattering (Waseda, 1980). The figure clearlyshows both EAM potentials are in close agreement with the tight binding andexperimental results. However, the main discrepancies are the peaks’ positionpredicted by both EAM those are slightly off. In addition S ( q ) of the mainpeak obtained by original EAM shows notably smaller value (12%).3.2 RelaxationEach cluster was relaxed as discussed in section 2.1 to minimize the poten-tial energy of the entire system at the beginning. For the small cub clusters, ize and shape-dependent melting mechanism of Pd nanoparticles 9 E ( e V / a t o m ) bcc ab initiofcc ab initiohcp ab initiobcc original EAMfcc original EAMhcp original EAMbcc optimized EAMfcc optimized EAMhcp optimized EAM Fig. 2
Validation of original and optimized EAM potentials for determining E c in compar-ison with ab initio calculation. a transformation to ico is expected during relaxation at room temperature(Baletto et al, 2002) The optimized EAM only predicts 2-cub transforma-tion to 2-ico while relaxation using original EAM shows 2-cub and 4-cub aretransforming to ico clusters in agreement with the previous result using thesame potential (Schebarchov and Hendy, 2006). This difference originates fromhigher energy barrier for the cub to ico transformation in optimized potential.As mentioned in the introduction the solid state transition is expected to occurvia an intermediate quasi-liquid state (Schebarchov and Hendy, 2006) whichseems to require higher energy using optimized EAM. In the following, we stickto ground state notations of the cluster and still call them 2-cub and 4-cub.We would like to remark that the cub structure is not the best fcc structurewith the lowest energy since its (100) facets are too large compared to (111)facets.The relaxation of higher shell number shows that there is no change in theshape and symmetry of the particle, indicating both potentials can successfullymodel the stable shapes of Pd clusters. S ( q ) Tight-binding [46]original EAMoptimized EAMExperiment [47]
Fig. 3
Validation of both EAM potentials for determining the S ( q ) of molten Pd at 1853 Kcompared with results of tight binding (Alemany et al, 1999) and experiment (Waseda, 1980)at same temperature. Fig. 4 illustrates the variation of U and C p versus T for the 8-cub particle(Pd ) and corresponding snapshots of point A – F for optimized EAM andG – L for original EAM. Both caloric curves present typical melting behaviorincluding an isotherm transition due to H mp . The gradual melting transition(non-isotherm melting) is more drastic in the curve obtained by original EAMwhich makes determining T mp harder. Nevertheless, T mp of 1277 and 1387 Kwere determined using main C p peaks for the original and optimized EAM, re-spectively. The difference came from the fact that optimized potential predictsbulk melting point more precisely and higher than original EAM. SnapshotsA – F belong to the optimized EAM which in general present more red atomswith higher potential energy. This is expected since U for the optimized EAMis always stand above original EAM values. It can be clearly seen in snapshotsA – D, the cub structure remains unchanged without any evidence of surfacemelting but slightly rounding in the corners. Snapshots E and F were takenat and after T mp those are indicating a homogeneous melting mechanism. Insnapshots G – I a solid-state cub to ico transition is evident (see the the movie ize and shape-dependent melting mechanism of Pd nanoparticles 11 in the supplementary materials). This is associated with a clear step changein U and corresponding local minimum in C p both at 1070 K from originalEAM. It has been shown previously that, a solid-state transition is accompa-nied with a step change in U and a local minimum in C p , which was referredto as negative C p (Zhang et al, 2010). The J – L snapshots were taken closeto melting point those showing a diagonal (corner to corner) melting whichstarted at the bottom and moving up where (111) facet on the top can still beseen. This is the reason for drastic non-isotherm melting in U of the originalEAM. Finally, snapshot L shows a complete liquid state.It is worth noting that the minor C p peak, is sometimes misinterpreted assurface melting (Pan et al, 2005). However, as mentioned in the introduction,the improper use of G ( r ) was the real reason for the misinterpretation. Nosurface melting is detected here as can be seen in snapshot H which is takenafter minor C p peak. Fig. 5 depicts the G ( r ) variation with temperature for 8-cub using both poten-tials. In both cases, the G ( r ) shows the same pattern for solid and liquid-likestates in agreement with previously reported patterns (Pan et al, 2005). Thesolid-like state at the bottom of each figure can be interpreted from 4 mainpeaks indicated by dashed line corresponding to 4 shells (not to be confusedwith 4th NNs). The 1st shell is consisting of 12 equidistant atoms and thuspresents a single peak while further shells contain more peaks due to theirgeometrical complexity. The liquid-like state at top of each figure consistingof 4 broad peaks those are also indicated by dashed lines. It can be clearlyseen that the increase in temperature causes peak broadening in both solidand liquid-like states. The melting occurs when there is step change in thepeaks’ positions (dashed lines). Fig. 5a shows homogeneous melting of all 4shells at about 1400 K using the optimized EAM in agreement with 1387 Kobtained from caloric curves. In Fig. 5b there are two step change in the peaks’positions at about 1000 and 1400 K. It is shown earlier that the first transitionbelongs to cub to ico transformation. However, since the step change in peaks’positions is more evident in the 4th peak (corresponding to 4th shell) it hasbeen misinterpreted as shell by shell melting in the previous study (Pan et al,2005). Here the optimized EAM is neglected since original EAM present a more com-plicated case with a solid-state transition and drastic non-isotherm melting.Fig. 6 presents the variation in the ratio of different structures with T accord-ing to the CNA for 8-cub cluster obtained by original EAM. The figure alsoincludes corresponding snapshots of the particle cross section at points A – Fthose indicated with dotted lines. The mirror change can be clearly seen infcc and disordered structure, including surface atoms, till ∼
940 K. The cross
A B CD E F400 600 800 1000 1200 1400 1600 1800 2000 2200 2400 T (K) -3.6-3.4-3.2 U ( e V / a t o m ) C p ( J / m o l . K ) A B CD E FG HI J K L U original EAM U optimized EAM Cp original EAM Cp optimized EAM G H IJ K L
Fig. 4
Variation of U and C p with temperature for 8-cub particle using optimized andoriginal EAM. (A – F) are the corresponding snapshots of points A – F obtained by optimizedEAM while (G – L) are snapshots of point G – L from original EAM. The color bar in theinset indicates potential energy level of each atom.ize and shape-dependent melting mechanism of Pd nanoparticles 13 Fig. 5
Variation of normalized G ( r ) with T for 8-cub cluster obtained by (a) optimizedand (b) original EAM force fields. The colorbar indicates normalized G ( r ) and dashed linesindicate its main peaks. section of the particle at point A in this region shows only nucleation of dis-ordered atoms close to the surface and remaining fcc atoms in green. In the940 – 1050 K range some hcp atoms are randomly appearing as can be seenin the snapshot B. It is worth noting that disordered phase exists everywhereeven in the core part before transition to ico. Then the transition occurs at1070 K, between point B and C, with a clear transition of fcc to disorderedatoms ( ∼ ∼ ∼
6% left disordered during the transition. Thus, CNA predicts an incompletetransition while we could not detect such thing using G ( r ) and caloric curves.Thereafter both fcc and hcp transform to disordered atoms as shown in theD – F snapshots. In specific, the diagonal melting is evident in D – F with adisordered lower half growing towards crystalline aggregates at the upper halfuntil total disorder is achieved.The local maximum in the ratio of disordered atoms, between B – C,was previously interpreted as “coexistence of solid-liquid phase” by focusingonly on the CNA criteria (Schebarchov and Hendy, 2006). However, our G ( r )
400 600 800 1000 1200 1400 T (K)1030507090 S t r u c t u r e r a ti o ( % ) A B C D E FfccdisorderedhcpA B CD E F
Fig. 6 variation of structure ratios with T and cross section of the corresponding snapshotsat A – F indicated in figure by dashed lines. result shows no liquid-like state during this transition. Thus the transitionoccurs with a so-called military diffusion mechanism, which is very fast andshort ranged, and the term “transition quasi-liquid state” is more preferredover solid-liquid coexistence. It is worth noting that the transition is limitedto 6 and 8-cub clusters possibly because for the larger clusters, the activationenergy required for this aim is higher than melting. ize and shape-dependent melting mechanism of Pd nanoparticles 15 Table 2
Dominant melting mechanism obtained using original (Foiles et al, 1986) and op-timized (Sheng et al, 2011) EAM force fields. The superscript (cid:63) denotes cub to ico transitionbefore melting and combined indicates simultaneous surface and diagonal melting mecha-nism. original EAM optimized EAM n cub ico cub ico2 uniform uniform uniform uniform4 uniform uniform uniform (cid:63) uniform6 diagonal (cid:63) diagonal surface diagonal8 diagonal (cid:63) diagonal surface diagonal10 surface diagonal surface diagonal12 surface diagonal surface diagonal14 surface diagonal surface diagonal16 surface diagonal surface combined18 surface diagonal surface combined20 surface diagonal surface combined n = 2 −
4) the ico is the only stable shape of clusters beforemelting those are melting uniformly due to the small diameter of particles. Forthe larger diameters, cub clusters are generally melting with surface meltingand ico clusters present diagonal melting. Optimized EAM predicts combinedsurface and diagonal melting for 16 – 20-ico clusters (10 −
12 nm in diameter).It is worth noting that there is also a difference in the surface melting of cuband ico clusters. The surface melting of the cub clusters is associated withnucleation of the liquid phase at (100) planes followed by its growth at thesurface which deforms polyhedral to a sphere with a lower surface area andsurface energy (see the related movie in the supplementary materials). Thisis interesting since it is shown experimentally that the (100) planes presentan incomplete surface melting in the bulk state while (110) and (111) planespresent complete and no surface melting, respectively (Vanselow and Howe,1988). Thus, the so-called two-stage melting is found to be associated withfaster growth of liquid phase at the surface before moving towards the core.While in the case of large ico clusters, 16 – 20-ico using optimized EAM, sur-face melting is in-fact simultaneous diagonal melting starting from differentcorners rather than formation of liquid shell followed by melting in the core(see the related movie in the supplementary materials). This difference can beexplained by more resistance of (111) planes against surface melting comparedto (100) ones as detected in many experiments. The interested reader is re-ferred to reference (Vanselow and Howe, 1988) and references therein. Thus,ico clusters with totally (111) planes at the surface present almost no surfacemelting. While cub cluster with both (100) and (111) planes at the surfaceshows surface melting. T mp T mb = 1 − (1 − q ¯ (cid:15) s (cid:15) v ) N s N t , q = ¯ Z s Z v (7)here N s stands for the number of surface atoms, q is the coordination numberratio with Z v equal to 12 for interior atoms and ¯ Z s as the average coordinationnumbers of surface atoms. The (cid:15) v and ¯ (cid:15) s respectively are bond energies ofinterior and surface atoms which the latter consists of cluster faces, edges andcorners.Fig. 7 compares normalized T mp obtained from caloric curves in comparisonwith previous MD (Pan et al, 2005; Miao et al, 2005) and MC (Lee et al, 2001)simulations. The figure also contains values calculated from Eq. 7 assuming thebond strength for the surface and interior atoms to be equal. It is evident that T mp values obtained by optimized EAM always stands between original EAMand Safaei model. In addition, previous MD results using SC force field arepredicting slightly higher values than original EAM but lower than optimizedEAM. The original EAM data set shows a clear change around 5 nm for the cubdue to the fact that 2 and 4-cub already transformed to ico during relaxation;and 4 – 5 nm cub (corresponding to 6 and 8-cub) are transforming to ico beforemelting. In the case of optimized EAM, T mp is higher for ico till 5 nm and islower for D > T mp for ico particles larger than 6 nmwhich means they are more stable than cub clusters. This might be associatedwith the γ sv difference in two potentials.3.6 Size-dependent melting enthalpyAttarian Shandiz and Safaei (2008) proposed the following model for the size-dependent H mp of clusters. H mp H mb = [1 − − q ) D D + D ] . [1 + 3 (cid:60) T mb H mb ln(1 − − q ) D D + D )] (8)with D being specific diameter which the entire atoms are located at thesurface (i.e. N s = N t ) and H mb is the bulk melting enthalpy. ize and shape-dependent melting mechanism of Pd nanoparticles 17 D (nm) T m p / T m b Fig. 7
Dependence of the normalized T mp on the size of particle obtained from U in com-parison with MD (Pan et al, 2005; Miao et al, 2005) and MC (Lee et al, 2001) simulation aswell as Safaei model (Safaei, 2010). The entire data sets were normalized to T mb = 1828 K(Vanselow and Howe, 1988). The variation of H mp with particle size is shown in Fig. 8 in comparisonwith Attarian Shandiz and Safaei (2008) model and previous MD results (Panet al, 2005; Miao et al, 2005). It is worth noting that, the model shows anegligible difference between ico and cub and thus it plotted for different ¯ Z s of6 and 3 corresponding to q = 0 . H mp with the particle size in agreementwith the result of the optimized EAM. While the original EAM underestimates H mp for particles larger than 5 nm. It can be seen SC potential (Pan et al,2005) fails to predict correct values for cub and ico by a huge error ( ∼ H mp < Z s of 6 and 3 meaning that melting smaller particlesare exothermic and favorable. This failure originated from the crystalline basisof the model, which is not defined for a few atoms. Both EAM potentialspredict a unified trend for ico particles. For the cub clusters, however, there isstep change at ∼ D (nm) H m ( k J / m o l ) cub original EAMico original EAMcub optimized EAMico optimized EAMEq. (8) q = 0.5Eq. (8) q = 0.25cub MD [18]ico MD [18]sphere MD [19]Bulk [52] Fig. 8
Variation of H mp with the diameter of Pd clusters obtained by optimized and originalEAM potential, in comparison with previous MD results (Pan et al, 2005; Miao et al, 2005)as well as Attarian Shandiz and Safaei (2008) model using D = 0 . In conclusion, the stability and melting behavior of Pd – Pd clusters ofcub and ico shape studied by molecular dynamics simulation using two EAMforce field i.e. original and highly optimized. Different melting criteria werediscussed those make a clear picture of the melting process. The result showssmall cub clusters (Pd – Pd ) are unstable those are melting uniformly dueto their small diameter. The only exception is Pd using optimized EAMwhich transforms to ico at elevated temperatures before melting. A similarcub to ico transition is predicted by original EAM for Pd – Pd clusterand for larger clusters both cub and ico are stable up to the melting point.In general, as detected by G ( r ) and CNA, the cub clusters present surfacemelting while ico clusters are melting diagonally (corner to corner) thanksto their partitioned structure. However, large ico cluster also present surfacemelting due to simultaneous diagonal melting from different corners. This isentirely different from surface melting observed in cub clusters. The later isassociated with nucleation of liquid phase at the (100) planes and its growthat the surface before moving inward. ize and shape-dependent melting mechanism of Pd nanoparticles 19 Conflict of InterestThe authors declare that they have no conflict of interest.
References
Alavi S, Thompson DL (2006) Molecular dynamics simulations of the melt-ing of aluminum nanoparticles. The Journal of Physical Chemistry A110(4):1518–1523Alemany MMG, Di´eguez O, Rey C, Gallego LJ (1999) Molecular-dynamicsstudy of the dynamic properties of fcc transition and simple metals in theliquid phase using the second-moment approximation to the tight-bindingmethod. Physical Review B: Condensed Matter 60(13):9208–9211, URL https://link.aps.org/doi/10.1103/PhysRevB.60.9208
Allen MP, Tildesley DJ (1989) Computer simulation of liquids. Oxford uni-versity press, OxfordAshcroft NW, Mermin ND (1976) Solid State Physics. Holt, Rinehart andWinston, New YorkAttarian Shandiz M, Safaei A (2008) Melting entropy and enthalpy of metallicnanoparticles. Materials Letters 62(24):3954–3956Baletto F, Mottet C, Ferrando R (2000) Reentrant morphology transition inthe growth of free silver nanoclusters. Physical review letters 84(24):5544Baletto F, Mottet C, Ferrando R (2001) Microscopic mechanisms of the growthof metastable silver icosahedra. Physical Review B 63(15):155,408Baletto F, Ferrando R, Fortunelli A, Montalenti F, Mottet C (2002) Crossoveramong structural motifs in transition and noble-metal clusters. The Journalof chemical physics 116(9):3856–3863Bertoldi DS, Mill´an EN, Guillermet AF (2017) Thermodynamics of the melt-ing process in au nano-clusters: Phenomenology, energy, entropy and quasi-chemical modeling. Journal of Physics and Chemistry of Solids 111:286–293Chen T, Zhang Y, Xu W (2016) Size-dependent catalytic kinetics and dynam-ics of Pd nanocubes: a single-particle study. Physical Chemistry ChemicalPhysics 18(32):22,494–22,502Dinsdale AT (1991) SGTE data for pure elements. Calphad 15(4):317–425Faken D, J´onsson H (1994) Systematic analysis of local atomic structurecombined with 3d computer graphics. Computational Materials Science2(2):279–286Foiles S, Baskes M, Daw M (1986) Embedded-atom-method functions for thefcc metals Cu, Ag, Au, Ni, Pd, Pt, and their alloys. Physical Review B33(12):7983Fu Q, Zhu J, Xue Y, Cui Z (2017) Size- and shape-dependent melting enthalpyand entropy of nanoparticles. Journal of Materials Science 52(4):1911–1918, DOI 10.1007/s10853-016-0480-9, URL https://doi.org/10.1007/s10853-016-0480-9
Goldstein AN, Echer CM, Alivisatos AP (1992) Melting in semiconductornanocrystals. Science 256(5062):1425–1427Iida T, Guthrie RIL (1988) The physical properties of liquid metals. ClarendonPress, OxfordJiang Q, Shi HX, Zhao M (1999) Melting thermodynamics of organic nanocrys-tals. The journal of chemical physics 111(5):2176–2180Jiang Q, Yang CC, Li JC (2002) Melting enthalpy depression of nanocrystals.Materials Letters 56(6):1019–1021Jos´e-Yacam´an M, Mar´ın-Almazo M, Ascencio JA (2001) High resolution TEMstudies on palladium nanoparticles. Journal of Molecular Catalysis A: Chem-ical 173(1):61–74Kateb M, Dehghani K (2012) Comparison of fracture behavior of sharp withblunt crack tip in nanocrystalline materials by molecular dynamics simula-tion. International Journal of Modern Physics: Conference Series 5:410–417Kraftmakher YA, Strelkov P (1970) Equilibrium concentration of vacanciesin metals. In: Seeger A, Schumacher D, Schilling W, Diehl J (eds) Vacan-cies and Interstitials in metals: International Conference Proceeding, NorthHolland, Amsterdam, p 59, held in J¨ulich, Germany, September 23-28, 1986Lee YJ, Lee EK, Kim S, Nieminen RM (2001) Effect of potential energy dis-tribution on the melting of clusters. Physical review letters 86(6):999Liang T, Zhou D, Wu Z, Shi P (2017) Size-dependent melting modes andbehaviors of Ag nanoparticles: a molecular dynamics study. Nanotechnology28(48):485,704Lindemann FA (1910) The calculation of molecular vibration frequencies.Physikalische Zeitschrift 11:609–612Miao L, Bhethanabotla VR, Joseph B (2005) Melting of Pd clusters andnanowires: a comparison study using molecular dynamics simulation. Phys-ical Review B 72(13):134,109Pan Y, Huang S, Liu Z, Wang W (2005) Molecular dynamics simulation ofshell-symmetric Pd nanoclusters. Molecular Simulation 31(14-15):1057–1061Plimpton S (1995) Fast parallel algorithms for short-range molecular dynam-ics. Journal of computational physics 117(1):1–19Plimpton SJ, Thompson AP (2012) Computational aspects of many-body po-tentials. MRS bulletin 37(5):513–521Poole Jr CP, Owens FJ (2003) Introduction to nanotechnology. John Wiley &Sons, New JerseyQi W (2016) Nanoscopic thermodynamics. Accounts of chemical research49(9):1587–1595Qi Y, C¸ aˇgin T, Johnson WL, Goddard III WA (2001) Melting and crystal-lization in Ni nanoclusters: The mesoscale regime. The journal of chemicalphysics 115(1):385–394Rangel E, Sansores E, Vallejo E, Hern´andez-Hern´andez A, L´opez-P´erez P(2016) Study of the interplay between N-graphene defects and small Pdclusters for enhanced hydrogen storage via a spill-over mechanism. PhysicalChemistry Chemical Physics 18(48):33,158–33,170 ize and shape-dependent melting mechanism of Pd nanoparticles 21
Rao CN, Rao KK (1964) Effect of temperature on the lattice parameters ofsome silver-palladium alloys. Canadian Journal of Physics 42(7):1336–1342Rossi G, Ferrando R (2007) Freezing of gold nanoclusters into poly-decahedralstructures. Nanotechnology 18(22):225,706Safaei A (2010) The effect of the averaged structural and energetic featureson the cohesive energy of nanocrystals. Journal of Nanoparticle Research12(3):759–776Safaei A, Shandiz MA, Sanjabi S, Barber ZH (2008) Modeling the meltingtemperature of nanoparticles by an analytical approach. The Journal ofPhysical Chemistry C 112(1):99–105, DOI 10.1021/jp0744681, URL http://dx.doi.org/10.1021/jp0744681
Schebarchov D, Hendy S (2006) Solid-liquid phase coexistence and structuraltransitions in palladium clusters. Physical Review B 73(12):121,402Schmidt M, Kusche R, von Issendorff B, Haberland H (1998) Irregularvariations in the melting point of size-selected atomic clusters. Nature393(6682):238–240Sheng HW, Kramer MJ, Cadien A, Fujita T, Chen MW (2011) Highly op-timized embedded-atom-method potentials for fourteen fcc metals. Physi-cal Review B 83(13):134,118, URL https://link.aps.org/doi/10.1103/PhysRevB.83.134118
Shim JH, Lee BJ, Cho YW (2002) Thermal stability of unsupported goldnanoparticle: a molecular dynamics study. Surface science 512(3):262–268Steinhardt PJ, Nelson DR, Ronchetti M (1983) Bond-orientational order inliquids and glasses. Physical Review B 28(2):784Stukowski A (2009) Visualization and analysis of atomistic simulation datawith OVITOthe open visualization tool. Modelling and Simulation in Ma-terials Science and Engineering 18(1):015,012Tsuzuki H, Branicio PS, Rino JP (2007) Structural characterization of de-formed crystals by analysis of common atomic neighborhood. Computerphysics communications 177(6):518–523Tyson WR, Miller WA (1977) Surface free energies of solid metals: Estimationfrom liquid surface tension measurements. Surface Science 62(1):267–276Vanselow R, Howe RF (1988) Chemistry and physics of solid surfaces VII,vol 10. Springer-Verlag Berlin Heidelberg, BerlinVerlet L (1967) Computer ”experiments” on classical fluids. I. Thermodynam-ical properties of Lennard-Jones molecules. Physical Review 159(1):98Wang W, Xu J, Zhao Y, Qi G, Wang Q, Wang C, Li J, Deng F (2017) Facetdependent pairwise addition of hydrogen over Pd nanocrystal catalysts re-vealed via NMR using para-hydrogen-induced polarization. Physical Chem-istry Chemical Physics 19(14):9349–9353Waseda Y (1980) The structure of non-crystalline materials: liquids and amor-phous solids. McGraw-Hill International Book Co., New YorkWestergren J, Nordholm S (2003) Melting of palladium clustersdensityof states determination by Monte Carlo simulation. Chemical physics290(2):189–209