Determination of melting mechanism of Pd24Pt14 nanoalloy by multiple histogram method via molecular dynamics simulations

Chemical Physics 388 (2011) 23–30

Contents lists available at ScienceDirect

Chemical Physics journal homepage: www.elsevier.com/locate/chemphys

Determination of melting mechanism of Pd24Pt14 nanoalloy by multiple histogram method via molecular dynamics simulations Hassan Yousefi Oderji a,b, Hongbin Ding a,⇑ a b

School of Physics and Optical Electronics, Dalian University of Technology, Dalian, Liaoning 116024, PR China School of Chemistry, University College of Science, University of Tehran, Tehran 14155, Iran

a r t i c l e

i n f o

Article history: Received 5 April 2010 In final form 8 July 2011 Available online 20 July 2011 Keywords: Pd24Pt14 nanoalloy Melting mechanism Multiple histogram Self time–space correlation function Caloric curve Heat capacity

a b s t r a c t Melting mechanism for the Pd24Pt14 nanoalloy has been determined by classical simulations in which quenching by the steepest descent method is coupled to the isothermal molecular dynamics and the nonergodicity of simulations is removed by the multiple histogram method. The Gupta many-body model is used for interatomic potentials. The melting characteristics are determined by the analysis of variations in the potential energy, the heat capacity, the vibrational density of states and the self time–space correlation functions with temperatures starting from 50 K up to 2200 K. The calculations indicate that the melting of Pd24Pt14 nanoalloy occurs at T = 747 K following structural transitions from the truncated-octahedral to the icosahedral basin of structures. A postmelting phenomenon has also been observed at temperatures between 900 and 1100 K, which is related to the homogeneous melting. This improves our knowledge about the phase transition features in 38-atom nanoalloys. Ó 2011 Elsevier B.V. All rights reserved.

1. Introduction Bimetallic nanoparticles or nanoalloys are new materials with interesting features in comparison with both pure nanometals and bulk alloys. They have optical, mechanical, and magnetic properties different from those of bulk alloys because of the size effect and from pure nanometals because of the synergistic effect. Moreover, these properties can be tuned by changing the composition, size, atomic order, and structure of particles. This especially makes them more powerful than the corresponding materials in their usage in the biomedical, catalytic, optical, and electronic industries [1–10]. Palladium and Platinum are of interest because they have been widely used as catalyst in the catalytic reduction of unsaturated hydrocarbons. Recently, Pd–Pt nanoalloys have been examined to increase the catalytic efficiency and to reduce the harmful side effects in hydrogenation reactions [1,11]. In a study of hydro-dearomatization catalysis by PdshellPtcore nanoalloys, Bazin et al. confirmed that the optimum Pd/Pt ratio is size-dependent and must be adjusted to maintain a Pt core surrounded by a complete Pd shell. Pd24Pt14 has a magic structure with high symmetry of the truncated octahedral (TO) and the optimum Pd/Pt ratio with a complete Pd shell [12]. The present work deals with the thermal stability of this nanoalloy. The significance of studying on thermal stability and the phase transition of Pd24Pt14, besides its importance to require the knowl⇑ Corresponding author. Tel./fax: +86 411 84706730. E-mail address: [email protected] (H. Ding). 0301-0104/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.chemphys.2011.07.011

edge for designing chemical reactions and exploring new applications, is due to the interesting features of phase transitions in 38-atom TO clusters. Studies on 38-atom Lennard-Joens (LJ) cluster have shown that unlike the most small LJ clusters which are based on an icosahedral core, LJ38 has a truncated octahedron structure as global minimum [13,14]. The icosahedral-core based structures, then, are at similar energies to on the potential surface, separated by a large energy barrier. This status affects both the simulation procedure and phenomenological properties of a system. Trapping the system in its global minimum basin or in the basin of the next highest energy isomer makes it a nonergodic system which needs to be handled with the appropriate simulation methods to break its nonergodicity [13,14]. For Pd24Pt14, the potential surface is even more complex than LJ38. According to the structural motifs reported by Paz-Borbón and coworkers, there are ten other symmetric structures with energies close to the global minimum [15]. Therefore, TO-structure of Pd24Pt14 is also expected to show especial melting features because of its particular shape of the potential surface. In order to study the melting process of Pd24Pt14, we used the steepest descent quenching method coupled to the isothermal molecular dynamics simulations to determine the mechanism of the phase transition, and we also analyzed the caloric curve and the heat capacity changes with temperature to find melting characteristics. More details about the dynamics of Pd24Pt14 nanoparticle in different temperatures have been obtained by using the velocity autocorrelation function (VACF) and the power spectrum of all constituent atoms. In addition, to have a superficial

24

H.Y. Oderji, H. Ding / Chemical Physics 388 (2011) 23–30

comparison between the vibrational frequencies of individual atoms and the frequencies of normal modes, we calculated the instantaneous normal mode (INM) intensities (including imaginary frequencies) by diagonalizing the dynamical matrix at some points in configuration space. We also used the self time–space correlation functions (SSTCF) to study the dynamics of structure evolution by calculating the values at a determined distance close to the each atom. In the field of phase transition, nanoalloys also have more complex behavior than the corresponding bulk alloys so that their phase transitions generally occur in a region of temperature rather than at a distinct sharp melting point and the dynamical phase coexistence phenomenon is one of the distinctive features of nanoparticles [16]. According to Berry and Smirnov [17], phase transitions in a cluster may be considered as being a result of configuration transformations. Hence, the phase transition mechanism of nanoparticles can be determined by exploring all configurations and their lifetime along the molecular dynamics (MD) phase space trajectories at solid, liquid and coexistence regions. However, the instant position of an atom in a trajectory is a resultant of its position in a frozen equilibrium configuration and the displacement caused by vibrations. We used the quenching steepest descent method to remove these displacements and to bring the atoms in their equilibrium positions. Coupling this method to MD simulation trajectories has been used to investigate the physical mechanism of transition from the rigid to nonrigid behavior of small Argon clusters by Beck and Berry [18]. Dynamics of the phase transition can be explored by studying the vibration density of state (VDOS) at different temperatures. Several methods are reported to calculate VDOS from the molecular dynamics simulation [19]. The Fourier transform of the velocity autocorrelation function [20], namely the power spectrum, and the instantaneous normal mode analysis [19,21] are among the most important methods. While Adams and Stratt [21] have introduced the instantaneous normal mode analysis as a probe of the cluster dynamics, VACF and the power spectrum have been reported as a method for calculating the phonon density of states and the vibrational properties of the metallic and nonmetallic clusters [22,23]. Calvo and Balbuena [23] calculated VDOS by the Fourier transforming of VACF for the different atomic species in Pt-Ag and Pt-Au clusters and explained the ability of VDOS measurements to obtain information about the chemical ordering in the experimentally produced nanoalloys. In the case of melting process, this method has been used by Jellinek and coworkers [24] to study the solid–liquid phase transition in the isoergic Ar13 cluster. Beck et al. [25] and Berry and Beck [26] developed this study to larger clusters Arn (n > 13). In all of these works, VACF is calculated by averaging over all cluster atoms. Recently, Hsu and coworkers [27] calculated VACF and its corresponding power spectrum for the individual atoms in a cluster and claimed that method is able to explore more details for locations and environments of individual atoms at different temperatures. In the other hand, SSTCF can reflect both the structural and the dynamical features of a cluster. SSTCF is a function that shows the correlation between the position r0 at time t0 and the position r at time t. It is expected to be a constant value for solids and a decaying value for liquids at the long correlation time when it is calculated at a close distance to each atom. The constant value of SSTCF at long correlation time is a result of anharmonic vibrations. Hence, SSTCF can be used as an indicator for the phase transition in nanoparticles. This quantity has been previously used by Battachara et al. [28] to study the dynamical phase coexistence in Ar13 cluster in a different approach from the present work. In this work, we investigate the ability of SSTCF as a dynamic criterion to elucidate melting process in nanoalloys in comparison with the caloric

curve, the constant volume heat capacity, and the vibrational density of states for Pd24Pt14. This paper is organized in four sections as follows: Section 1 presents an introduction to nanoalloys and the phase transition in nanoparticles. Simulation details, the Gupta potential model and the multiple histogram method are explained in Section 2. Section 3 deals with the analysis and discussion of melting mechanism using the steepest descent quenching method, the caloric curve, the heat capacity, the vibrational density of states and the self time–space correlation functions. Finally, the results of all criteria are compared and concluded in Section 4.

2. Simulation details Recently, the Gupta empirical-potential description of Pd-Pt clusters has been found to be reasonably accurate by Paz-Borbón et al. [15], via comparing the results of this model from the genetic algorithm with the results from the density functional theory (DFT). By this method, they found 11 stable configurations for Pd24Pt14 in which TO structure has lowest configuration energy, and it is known as global minimum (GM) of Pd24Pt14. In this structure, Pt atoms construct the octahedral cluster core with a perfect growth of Pd atoms on top of the (111) faces so that the particle surface presents eight hexagonal and six square faces, and the grown cluster holds its Oh symmetry. The octahedral core and the whole nanoparticle are shown in Fig. 1(a) and (b). Starting from TO structure at low temperatures, the molecular dynamics simulations of this work were carried out using the DLPOLY molecular dynamics parallel simulation package [29] in three stages. In all of the stages, the molecular dynamics simulations were performed by the velocity–Verlet algorithm with a time step of 1 fs in NVT ensemble with the Evans thermostat [30]. In short time simulation stage, we started with initial temperature of 50 K and increased the temperature of the system by rate of 1 K/ns up to the cluster evaporation which occurs beyond T = 2200 K. This stage was repeated three times to produce three configurations at every temperature as the required input data for next stages. In the second stage, we performed long time simulations at some specific temperatures for increasing the statistical precision and sampling from all of the phase space, especially in the coexistence phase region. We examined the effect of the simulation length on the results by values of 25, 50, 100, and 200 ns. The instant velocities and positions of all particles were sampled every 5 fs in the second and third stages. The simulations of the second stage were performed one time for temperatures between 300 K and 2200 K and were repeated two more times for temperatures corresponding to the phase coexistence region. Nevertheless, as it will be explained later, the extension of the simulation length, even up to 200 ns, cannot improve the statistical precision of the results and the heat capacity values still show strong fluctuations in the coexistence phase region. Furthermore, the mechanism of phase

Fig. 1. The octahedral Pt14 core (a) and the truncated octahedron (TO) structure of Pd24Pt14 nanoalloy (b).

25

H.Y. Oderji, H. Ding / Chemical Physics 388 (2011) 23–30

2.1. The Gupta potential model

N X N N X 1X V ij ðr ij Þ þ Fðqi Þ; 2 i¼1 j–i i¼1

ð1Þ

where F(qi) represents the energy of embedding an atom in the bulk density, qi, which is defined as:

qi ¼

N X

qij ðrij Þ;

qij ðrij Þ ¼

n2a;b

exp 2qa;b

pffiffiffiffiffi Fðqij Þ ¼  qi :

r ij  r 0a;b

;

ð4Þ ð5Þ

In these equations, rij represents the distance between atoms i and j; a and b subscripts refer to the element type. The parameters of A, p, q, r0, and n used in this study are from Paz-Borbón et al. [15]. 2.2. Multiple histogram method In a range of temperatures, the Pd24Pt14 nanoalloy with TO structure is not an ergodic system, so its thermodynamic parameters, obtained by the classical simulations, strongly fluctuate and should be determined by appropriate techniques to breakthrough nonergodicity. Multiple histogram method is a technique to calculate thermodynamic properties of a system at desired temperatures (or pressures, etc.) using the data from Monte Carlo or molecular dynamics simulations at other temperatures (or pressures, etc.) [31,32]. This method is useful especially when a system is not ergodic at some temperatures. In the multiple histogram method, the unnormalized probability of the energy E at the inverse temperature of b = 1/(kBT) can be estimated from the probability of energy E at other temperatures as follows [32]:

PR

PðE; bÞ ¼

1 n¼1 g n N n ðEÞ expðbEÞ ; PR 1 m¼1 g m N m;tot expðbm E  fm Þ

2

4

8

6

10

12

E, (eV) Fig. 2. The normalized probability of P(E) at three different temperatures related to the TO basin (T = 590 K), the icosahedral basin (T = 1100 K), and the phase coexistence region (T = 750 K).

tion; fm is the free energy of system at bm, which is related to unnormalized probability at bm as follows:

expðfm Þ ¼

X

PðE; bm Þ:

ð7Þ

E

hAi ¼

! ;

0

ð3Þ

!

r 0a;b r0a;b

0.01

ð2Þ

and Vij(rij), qij(rij), and F(qi) are defined as:

r ij  r 0a;b

0.015

We used the direct iteration method to solve the Eq. (6) at all simulation temperatures to find the values of fm parameters as in Eq. (7). Once the probabilities of P(E, b) are known, the average values of any function of configuration space, namely A, can be obtained at inverse temperature of b as follows:

j¼1;j–i

V ij ðrij Þ ¼ Aa;b exp pa;b

0.02

0.005

In the molecular dynamics simulations of this work, the Gupta potential model has been used to calculate the forces on atoms. This model is a semi-empirical many-body potential model based on the second moment approximation to tight-binding theory [33] which describes the bonding energy of a metal atom in terms of the local electronic density. In this model, the total potential energy, E, is given by the sum over the both bonding and repulsive energies of all atoms as follows:

E¼

T=590 K T=750 K T=1100 K

0.025

Normalized P(E)

transition obtained by the steepest descent quenching method shows that the system is not ergodic in temperature range between 600 and 900 K. Therefore, we used the multiple histogram method [31,32] to overcome the nonergodicity. In order to that, we simulated the system using the output data of the first stage in the temperature range between 50 and 2200 K at every 10 K for 10 ns and estimated the free energies of the system at these temperatures using direct-iteration solution of equations which are presented in Section 2.2

ð6Þ

where gn = 1 for the independent configurations and gn = 1 + 2sn, otherwise; sn is the correlation time; R is the number of computer simulations, and Nn(E)/Nn,tot is the ratio of configurations with energy E to all simulated configurations at bn in nth computer simula-

X i

PðEi ; bÞ AðiÞ P : PðEj ; bÞ

ð8Þ

j

Fig. 2 shows the normalized values for P(E) at temperatures of 590, 750, and 1100 K. At temperatures lower than 590 K, the system is in TO basin of energies, but at higher temperatures, it can jump to the basin of icosahedral structures. At T = 750 K, the probability of finding system in these two basins is equal, and finally at around T = 1100 K, the probability of finding a configuration belong to the TO basin is almost zero. Hence, the Pd24Pt14 nanoalloy with TO structure is not an ergodic system in temperatures between 590 and 1100 K. 3. Mechanism of phase transition After simulation, phase trajectories should be analyzed to obtain the required information about the melting temperature and mechanism of the nanoparticle. We used the steepest quenching method to determine the most stable configurations at different temperatures, and the variations in the caloric curve and the heat capacity to determine the melting temperature. We also used VDOS and SSTCF dynamical criteria to insure the diffusive motions of all individual atoms in icosahedral basin of structures. 3.1. Steepest descent quenching method In order to understand the mechanism of transition from the solid-like to the liquid-like phases of Pd24Pt14, we used the steepest descent energy minimization method to quench instantaneous atomic configurations from the molecular dynamics simulations. A simulated instantaneous configuration is a snapshot of oscillating atoms; therefore, we need to quench the oscillations completely to find the nearest local potential minimum in which the system oscillates around it. The local potential minima on the energy landscape,

26

H.Y. Oderji, H. Ding / Chemical Physics 388 (2011) 23–30

Fig. 4. Excited configuration of Pd24Pt14 nanoalloy in coexistence phase region; the icosahedron Pt atoms of core (a), and the whole configuration with pentagonal faces (b). The values of configuration energies are relative to the global minimum.

Fig. 3. Some excited configurations of Pd24Pt14 nanoalloy simulated at temperatures lower than T = 700 K. The values of configuration energies are relative to the global minimum.

which are obtained by this method, and their transition pathways can give a mechanistic description of the phase transition. In this method, the following differential equation,

r_ ¼ rE;

ð9Þ

Relative configuration energy, ΔEc (eV)

is solved for some selected configurations on MD simulated trajectory at determined time intervals [18]. In order to access all local minima in a trajectory, the time interval should be adjusted with the temperature. At low temperatures, the time interval can be set to several times of vibrational periods while it should be shorter than one period at higher temperatures i.e., in the coexistence and the liquid-like regions. In our simulations, quenches were performed every 5 ps for temperatures lower than 600 K. All of these quenches lead to the

same TO structure (Fig. 1(b)), indicating to stability of the global minimum up to this temperature. The first excited configuration occurred at around T = 600 K at which one of the hexagonal faces slightly slides along a line parallel to its vicinal square faces (Fig. 3(b)). The lifetime of this configuration was less than 1 ps. Some of the more excited short lifetime configurations were also observed at higher temperatures (Fig. 3(c)–(d)), but the first stable excited configuration occurred when the octahedral structure of core transformed to icosahedral and the hexagonal arrangement of some shell atoms changes to the pentagonal. The most symmetric structure of the icosahedral basin of structures is presented in Fig. 4. The relative configuration energy of Pd24Pt14 to the energy of global minimum is presented in Fig. 5 as a function of time for the trajectories at different temperatures after quenching by the steepest descent method. Fig. 5(a) shows that the particle is in its global minimum for the most of simulation time (100 ns) at T = 695 K. Three different simulations at T = 700 K (Fig. 5(b–d)) indicate that the system is not ergodic at a region of temperatures. At this temperature system is mostly in equilibrium between the TO global minimum basin and the icosahedral basin of structures. However, Fig. 5(c) at time between 50 and 60 ns shows a TO structure for the particle in which one of the Pt core atoms is displaced with one of the Pd shell atoms. This and other possible kinds of

3.0 2.0 1.0 0.0

(a)

3.0 2.0 1.0 0.0

(b)

3.0 2.0 1.0 0.0

(c)

3.0 2.0 1.0 0.0

(d)

3.0 2.0 1.0 0.0

(e)

0

20

40

60

80

100

Time, t (ns) Fig. 5. Relative configuration energy of Pd24Pt14 nanoalloy obtained by Steepest descent quenching of instant configurations along molecular dynamics simulated trajectories at T = 695 K (a), T = 700 K (b, c, d), and T = 1000 K (e).

27

H.Y. Oderji, H. Ding / Chemical Physics 388 (2011) 23–30

displacements lead to the contribution of all atoms in the liquidlike behavior of system at this stage of melting. Fluctuations in the thermal properties of the nanoparticle in this region of temperatures are due to random nature of finding the cluster in the solid or liquid phases during limited-time simulations. 3.2. Caloric curve The caloric curve and the heat capacity are two common thermodynamical criteria which are often used to identify the melting transition [1,16,34–36]. The caloric curve reports the average energy of the cluster as a function of temperature. For the completely solid or liquid phases, the caloric curve is a quite gradual increasing function of temperature while for coexistence region it usually presents a drastically change that is related to the latent heat. Fig. 6 presents the caloric curve of Pd24Pt14 nanoalloy. A sharp jump is observed at approximately T = 750 K. Furthermore, a weak trace of postmelting phenomenon can be deduced from the caloric curve at temperatures approximately between 850–1000 K. We also extended the length of simulations, particularly in the coexistence phase region, to check if the samplings are from all of the configuration space. Cross symbols in Fig. 6 have been obtained from 50 ns simulation length while Circles have been obtained from 25 ns. Although, the results are compatible in both low and high temperatures, but they have strong fluctuations at a rage of middle temperatures. Extending the length of simulation up to 200 ns was not efficient to smooth the results of this region, so these results have not been shown in Fig. 6. These observations are consistent with the phase transition mechanism and confirm that the system is nonergodic. The solid line in Fig. 6 has been obtained from the multiple histogram method to cover this region. The latent heat of melting can be estimated from the gap in energy at melting temperature, indicated by q in Fig. 6. More details about the melting can be obtained from further analysis using the heat capacity and the dynamical criteria.

C v ðTÞ ¼

3 hE2 i  hE2 i NkB þ ; 2 kB T 2

ð10Þ

where E is the potential energy, N is the number of atoms, and kB is the Boltzmann factor. In contrast to bulk systems, which generally exhibit singularity in the heat capacity versus temperature curves at the melting point, nanoparticles show one or more softened peaks in a broader region of temperature [16]. Fig. 7 shows the change of heat capacity values with temperature for the Pd24Pt14 nanoalloy, which are obtained from the different simulation lengths of 25 ns (Circles), 50 ns (Crosses), and 200 ns (Diamonds). At low and high temperatures, the results from both 25 ns and 50 ns are the same, but at temperatures between 700–1000 K, even simulations with 100 ns of duration show strong fluctuations. These fluctuations are due to the nonergodicity and the extension of the length of simulations to 200 ns was still inadequate to reduce them. The results from 100 ns are excluded from Fig. 7 because their uncertainties are similar to the results from 200 ns of simulations. However, the heat capacity values show a drastic change, at about T = 750 K, that can be related to the phase transition, the exact determination of this peak and subsequently the melting temperature needs to use an appropriate technique to breakthrough the nonergodicity for smoothing the results. The solid line in Fig. 7, obtained using the multiple histogram method, shows that the melting occurs at T = 747 K. Fig. 2 demonstrate that the system at 747 K is in equilibrium between two energy basins on potential surface. The basin with lower energy is related to the TO structure and the higher one is related to the icosahedral structures. Fig. 7 also shows a postmelting stage, at temperatures between 890 and 1100 K. At this stage, the system has energy values higher than the icosahedral basin and the homogeneous melting occurs. These results are in agreement with those of caloric curve and the phase transition mechanism. 3.4. Vibrational density of states

3.3. Heat capacity The constant volume heat capacity, Cv, is another thermodynamical criterion to study the melting process in nanoparticles. In canonical ensemble, this quantity is related to the potential energy fluctuations as follows [37]:

Vibrational density of states can be calculated from the molecular dynamics simulations by several methods, e.g., from the equations of motion with the forces given in the harmonic approximation [38], from the diagonalization of the dynamical matrix in the instantaneous normal mode analysis [19] and from the Fourier transform of velocity autocorrelation functions of all atoms, i.e., the power spectrum as follows [24]:

20

8 7 6

Cv (kB)

ΔEconf (eV)

15

10

5

5 4

q

3

0 0

500

1000

1500

2000

T (K) Fig. 6. Caloric curve for the Pd24Pt14 nanoalloy. Cross symbols have been obtained from 50 ns simulation length, and Circles are from 25 ns. The solid line is from the multiple histogram method and dash lines are extrapolation of the frozen and molten branches, and q is related to the latent heat.

0

500

1000

1500

2000

T (K) Fig. 7. Heat capacity values of Pd24Pt14 nanoalloy at different temperatures. Cross symbols have been obtained from 50 ns simulation length; and Circles are from 25 ns, and Diamonds are from 200 ns. The solid line is from the multiple histogram method.

28

H.Y. Oderji, H. Ding / Chemical Physics 388 (2011) 23–30

XðiÞ ðxÞ ¼ 2

Z

1

C ðiÞ ðtÞ cosðxtÞdt;

ð11Þ

0

where C(i)(t) is the normalized velocity autocorrelation function of ith atom [39]:

C ðiÞ ðtÞ ¼

hv i ðt0 Þ:v i ðt þ t 0 Þi ; hv i ðt 0 Þ:v i ðt0 Þi

ð12Þ

1 @2E Dijab ¼ pffiffiffiffiffiffiffiffiffiffiffiffi M i M j @r ai r bj

! ð13Þ

; fr0i g

for the configurations at every 5 ps snapshots of each simulated trajectory and then by diagonalizing dynamical matrices to obtain their eigenvalues, ks, where s stands for the vibration mode, and eigenvalues are related to vibrational frequencies by ks = (2pts)2. The symbol E in Eq. 13 is the potential energy; M is mass of the particle; i and j are atoms indices; a and b stand for x, y and z coordinates, and fr0i g denotes the current configuration [19]. The INM analysis is presented in Fig. 9 at some temperatures, near to the average temperatures as in Fig. 8, to have a superficial comparison between the

0.25 0.20 0.15 0.10 0.05 0.00

=407 K <ΔE>=1.99 eV

=711 K <ΔE>=3.66 eV

=728 K <ΔE>=3.75 eV

0.15

=682 K <ΔE>=4.68 eV

=723 K <ΔE>=4.82 eV

=805 K <ΔE>=5.47 eV

=915 K <ΔE>=6.50 eV

=985 K <ΔE>=7.33 eV

=2322 K <ΔE>=19.3 eV

(i)

2

Ω (ω) [Å /ps]

where vi is the velocity of ith atom, and the angular brackets indicate ensemble averaging over short sub-trajectories of the system with time origins of t0. Power spectra represent the underlying frequencies of the individual atoms [27]. The low frequency limit of the power spectrum is related to diffusive motions. Hence, the increase in the intensity of low frequencies indicates to the increase of diffusive motions. This is one of the characteristics for fluids, so it can be used to distinguish between solid and liquid phases. Velocities of atoms have to be rescaled artificially during the isothermal simulations to conserve temperature. To avoid the ambiguities in velocity autocorrelation functions due to rescaling of velocities, we performed the microcanonical simulations at some selected energies to calculate the values of power spectra, which are presented in Fig. 8. At the energies lower than 4.68 eV, the 38 curves of the individual atoms stay in 3 sets. This indicates that there are three kinds of atom-locations in TO structure of Pd24Pt14 nanoalloy, i.e. two different locations for Pt atoms (red lines) and one location for Pd atoms. TO structure withstands energy increasing up to 4.68 eV in which one of the Pt atoms shows a different power spectrum than others. The positions of the corresponding peaks for this atom can be compared with those for the central Pt atom in Pd20Pt13, which we found them at 234, 136 and 67 cm1 by simulation of Pd20Pt13 nanoalloy at T = 700 K. There are three kinds of locations for Pt atoms in Pd24Pt14 based on the icosahedral structures: central Pt atom, 12 vertices Pt atoms and the one on the surface of Pt13 icosahedron, but Fig. 8 at DE = 4.68 eV shows only two kinds of power spectra, indicating to only two locations. This is due to frequently replacements of the Pt atom on surface with the Pt atoms at vertices, which pro-

vides the same conditions for them in average. Thus, atoms have diffusive motions and system behaves as liquid while the overall icosahedral symmetry of structures is conserved. Furthermore, Fig. 8 shows two and five curves respectively at DE = 4.82 and DE = 5.47 eV with characteristic frequencies of the central Pt atom. This is due to replacement of one Pt atom at DE = 4.82 eV and four Pt atoms at DE = 5.47 eV with the central Pt atom during simulations. Therefore, it can be deduced that at sufficient time, all atoms would contribute to the diffusive motions and liquid behavior of the particle. At energies beyond DE = 7.33 eV, the intensity of the zero-limit frequencies increases with energy, which indicates to the increase of diffusive motions. However, it should be recalled that the system is nonergodic at energies of 4.68, 4.82 and 5.47 eV in Fig. 8 and the presented results are from the limitedtime simulations. The power spectra of VACF only produce the positive frequencies of the individual atoms while imaginary frequencies of VDOS may be also interesting. These frequencies are associated to the cross of barriers or saddle points in isomerizations; therefore, the ratio of these modes to the total modes is an increasing function with temperature. The drastically changes in this function occurs at the phase transitions. Adams and Stratt used this method as a probe of cluster dynamics [21]. Here, we performed instantaneous normal mode analysis by calculating the 3N  3N(N = 38) elements of the dynamical matrix,

0.10 0.05 0.00 0.10

0.05

0.00

0

100 200 300 400

0

100 200 300 400

0

100 200 300 400

-1

ω (cm ) Fig. 8. Power spectra of all constitute atoms of Pd24Pt14 nanoalloy in different energies. Red lines relate to Pt, and blue ones relate to Pd atoms. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

29

H.Y. Oderji, H. Ding / Chemical Physics 388 (2011) 23–30

T=695 K

T=400 K

T=700 K (c)

0.005

Arbitrary units

0 T=700 K (d)

T=750 K

T=800 K

T=900 K

T=1000 K

T=2200 K

0.005

0

0.005

0

-100 0

200

400

-100 0

200

400

-100 0

200

400

-1

Frequency (cm ) Fig. 9. The instantaneous normal mode analysis of Pd24Pt14 at different temperatures.

frequencies from VACF and frequencies from INM analysis. Fig. 9 shows that there are some imaginary frequencies, shown as negative values, even at low temperatures; these frequencies are due to the anharmonic vibrations [21]. However, during and after melting, the imaginary frequencies are related to the transition states in isomerization reactions pathways. Although, it is possible to exclude the anharmonicity contribution in the imaginary frequencies by using the quenched instead of the original configurations, but it has not been performed in the present work. Applying of this method needs the time interval to be very short to capture adequate number of transition states to have a reasonable statistical precision.

3.5. Self time–space correlation function Self time–space correlation (or Van Hove) function, Gs(r, t), measures the probability that an atom is at position r at time t given that the same atom was at the origin r = 0 at the initial time t = 0 [39]. This function can be evaluated from the simulation data for each atom as follows:

GsðiÞ ðr; tÞ ¼

M X 1 d½r  jr i ðtÞkÞ  r i ðt k þ tÞj; 4pr 2 M k¼1

ð14Þ

0.1 0.08 0.06 0.04

T=400 K

T=695 K

T=700 K (c)

T=700 K (d)

T=750 K

T=800 K

T=900 K

T=1000 K

T=2200 K

0.02

i

0.06

2

4πr Gs (r,t)

0 0.08 0.04 0.02 0 0.08 0.06 0.04 0.02 0

0

1

2

3

4

5

0

1

2

3

4

5

0

1

2

3

4

5

Time (ps) Fig. 10. Self time–space correlation functions of all constitute atoms of Pd24Pt14 nanoalloy at different temperatures. Red lines relate to Pt, and blue ones relate to Pd atoms; ‘c’ and ‘d’ labels refer to labels in Fig. 5. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

30

H.Y. Oderji, H. Ding / Chemical Physics 388 (2011) 23–30

where index k runs over a total of M selected time origins from a system at equilibrium, and d is the Dirac delta symbol. In order to study the thermal behavior of Pd24Pt14 nanoalloy, we calculated the GðiÞ s ðr; tÞ values at a constant distance (r = 0.2 Å) around all atoms. This is especially useful to obtain information about structural details of a nanoparticle as well as the thermal stability. The results of 1 atoms (red color curves) and 24 Pd atoms (blue color curves) at the same temperatures as in Fig. 9 are presented in Fig. 10. This figure shows that they coincide with each other in three sets at temperatures lower than T = 700 K and in two sets approximately at temperatures higher than T = 900 K. At temperatures between T = 700 K and T = 900 K, SSTCF behaves somewhat different for the same atoms. At origin time, there is not any autocorrelation between the current position of an atom and positions on the surface of its surrounding sphere, with radius r = 0.2 Å. A short while later, atoms will pass through the surface of spheres and correlations will be maximum. For diffusing motions, the probability that one atom would return to the surface is low, so SSTCF decays to zero. The rate of decays indicate to the ratio of diffusive motions. In contrast, for vibrating motions, one atom returns on the surface after a period of vibration. For Fig. 10 at T = 695 K, the constant values for SSTCFs at long time are due to anharmonicity. For temperatures lower than 700 K, all SSTCF curves tend toward or fluctuate around constant values. This indicate to the vibrational motions of atoms around equilibrium positions and then system is in a solid phase. Since liquids are amorphous and the same atoms would be expected to show the same behavior, this suggests that the homogeneous melting occurs at temperature beyond T = 900 K. At temperature between 700 and 900 K, SSTCF curves decay slowly, so all atoms have diffusive motions. However, the rate of these motions are low and their decays occur at long times. The slight differences in behavior of similar atoms in this region are due to the nonergodic nature of system. At thermodynamic scale of time, they would behave similarly in average. Fast decays occur at the high temperatures at which the probability of the backward motions is low and diffusivity is high. Thus, SSTCF can be used as a criterion to identify the phase status and the phase transitions in nanoparticles. 4. Conclusion The melting mechanism of Pd24Pt14 nanoalloy has been studied by the isothermal molecular dynamics simulations. The simulated configurations were observed to be as jumping between TO-structure basin and icosahedral basin of structures at some temperatures. Moreover, strong fluctuations in caloric and heat capacity curves are the evidences for the nonergodic behavior of system at these temperatures. Therefore, we used the multiple histogram method to overcome the nonergodicity and found the melting occurs at T = 747 K. Analysis of power spectra and SSTCF indicates that the transition from TO-structure basin of energies to the icosahedral basin of structures is not a solid–solid transition, and all atoms in icosahedral basin of structures present diffusive motions and contribute to the melting process. The peak in heat capacity curve with temperature show a postmelting shoulder between T = 900 and 1100 K. This is related to the homogeneous melting of the nanoalloy. In addition, a comparison between the vibrational frequencies of individual atoms, obtained by the Fourier transform of velocity

autocorrelation functions, and the frequencies of normal modes, obtained by diagonalizing dynamical matrices on configuration space, indicates that the power spectra are more useful than the instantaneous normal mode analysis method to obtain the details about the structural evolution of Pd24Pt14 while INM results with the imaginary frequencies are related to the transition states in isomerization pathways. Acknowledgments This work was supported by the Scientific and Technical Foundation of Liaoning Province (No. 20082168), the National Science Foundation of China (No. 10875023), Scientific and Technical Key Project (No. 108034) and Ph.D research program (No. 200801411040) of the Educational Ministry and National Magnetic Confinement Fusion Science Program of China (No. 2009GB106004, 2008CB717801). The authors gratefully acknowledge the support from the Research Council of Tehran University and the Plasma Physics Computational Center of Dalian University of Technology. H.Y acknowledges the Chinese Scholarship Council and the Iran Nanotechnology Initiative Council from their supports. References [1] R. Ferrando, J. Jellinek, R.L. Johnston, Chem. Rev. 108 (2008) 845. [2] U. Kreibig, M. Vollmer, Optical Properties of Metal Clusters, Springer Series in Material Science, vol. 25, Springer, Berlin, 1995. [3] E. Cottancin, J. Lermé, M. Gaudry, M. Pellarin, J.-L. Vialle, M. Broyer, B. Prével, M. Treilleux, P. Mélinon, Phys. Rev. B 62 (2000) 5179. [4] M. Haruta, Catal. Today 36 (1997) 153. [5] S. Lee, C. Fan, T. Wu, S.L. Anderson, J. Chem. Phys. 123 (2005) 124710. [6] A.M. Molenbroek, S. Haukka, B.S. Clausen, J. Phys. Chem. B 102 (1998) 10680. [7] L.M. Molina, B. Hammer, Phys. Rev. Lett. 90 (2003) 206102. [8] H. Portales, L. Saviot, E. Duval, M. Gaudry, E. Cottancin, M. Pellarin, J. Lermé, M. Broyer, Phys. Rev. B 65 (2002) 165422. [9] G. Rossi, A. Rapallo, C. Mottet, A. Fortunelli, F. Baletto, R. Ferrando, Phys. Rev. Lett. 93 (2004) 105503. [10] M. Valden, X. Lai, D.W. Goodman, science 281 (1998) 1647. [11] B. Coq, F. Figueras, J. Mol. Catal. A 173 (2001) 117. [12] D. Bazin, D. Guillaume, Ch. Pichon, D. Uzio, S. Lopez, Oil Gas. Sci. Tech. Rev. IFP 60 (2005) 801. [13] F. Calvo, J.P. Neirotti, D.L. Freeman, J.D. Doll, J. Chem. Phys. 112 (2000) 10350. [14] J.P. Neirotti, F. Calvo, D.L. Freeman, J.D. Doll, J. Chem. Phys. 112 (2000) 10340. [15] L.O. Paz-Borbón, R.L. Johnston, G. Barcaro, A. Fortunelli, J. Chem. Phys. 128 (2008) 134517. [16] A. Proykova, R.S. Berry, J. Phys. B-At. Mol. Opt. 39 (2006) R167. [17] R.S. Berry, B.M. Smirnov, Low Temp. Phys. 35 (2009) 256. [18] T.L. Beck, R.S. Berry, J. Chem. Phys. 88 (1988) 3910. [19] D. Caprion, H.R. Schober, J. Chem. Phys. 114 (2001) 3236. [20] J.H. Dickey, A. Paskin, Phys. Rev. 188 (1969) 1407. [21] J.H. Adams, R.M. Stratt, J. Chem. Phys. 93 (1990) 1332. [22] A. Kara, T.S. Rahman, Phys. Rev. Lett. 81 (1998) 1453. [23] S.R. Calvo, P.B. Balbuena, Surf. Sci. 581 (2005) 213. [24] J. Jellinek, T.L. Beck, R.S. Berry, J. Chem. Phys. 84 (1986) 2783. [25] T.L. Beck, J. Jellinek, R.S. Berry, J. Chem. Phys. 87 (1987) 545. [26] R.S. Berry, T.L. Beck, H.L. Davis, J. Jellinek, Solid–Liquid Phase Behavior in Microstructures, John Wiley & Sons, Inc., 2007. [27] P.J. Hsu, J.S. Luo, S.K. Lai, J.F. Wax, J.-L. Bretonnet, J. Chem. Phys. 129 (2008) 194302. [28] A. Bhattacharya, B. Chen, S.D. Mahanti, Phys. Rev. E 53 (1996) R33. [29] W. Smith, T.R. Forester, J. Mol. Graphics 14 (1996) 136. [30] D.J. Evans, G.P. Morriss, Comput. Phys. Rep. 1 (1984) 297. [31] A.M. Ferrenberg, R.H. Swendsen, Phys. Rev. Lett. 63 (1989) 1195. [32] T. Bereau, R.H. Swendsen, J. Comput. Phys. 228 (2009) 6119. [33] F. Cleri, V. Rosato, Phys. Rev. B 48 (1993) 22. [34] Z. Kuntová, G. Rossi, R. Ferrando, Phys. Rev. B 77 (2008) 205431. [35] A. Lyalin, A. Hussien, A.V. Solov’yov, W. Greiner, Phys. Rev. B 79 (2009) 165403. [36] A. Aguado, M.F. Jarrold, Annu. Rev. Phys. Chem. 62 (2011) 151. [37] D. Sabo, C. Predescu, J.D. Doll, D.L. Freeman, J. Chem. Phys. 121 (2004) 856. [38] C. Oligschleger, J.C. Schon, J. Phys.: Cond. Matt. 9 (1997) 1049. [39] J.M. Haile, Molecular Dynamics Simulation: Elementary Methods, 1st ed., John Wiley & Sons, Inc., New York, 1992.