Melting of palladium clusters—Canonical and microcanonical Monte Carlo simulation

Abstract

We present Monte Carlo simulations of single palladium clusters of 13, 34, 54, 55, 147 and 309 atoms. The clusters are modeled by a many-body potential and they have been simulated at constant temperature or constant total energy. The caloric curves of the clusters, with the exception of Pd₃₄, exhibit an S-bend at melting which is typical for a finite system. We have also observed the typical coexistence region of solid and molten clusters both in the canonical and the microcanonical ensembles. Pd₃₄, in contrast, melts without an accompanying peak in heat capacity and at melting the atoms become mobile without any significant change in geometric structure. For the larger clusters a free energy barrier inhibits phase switching. In some cases of phase change from molten to solid structure the barrier is of purely entropic character. By a conversion of the results in the microcanonical simulations into temperature-dependent data, the simulations at fixed temperature and fixed total energy have been compared. The agreement is in most cases good. The results are furthermore compared to earlier molecular dynamics simulations with the Nosé–Hoover thermostat. These results are in good agreement with the Monte Carlo simulations as well.

---

I. Introduction

During the last three decades, keen interest has been shown in investigation of the melting of clusters, mainly computationally but also experimentally. The reasons are manifold. In the field of experimental cluster research, a knowledge of the thermal state of free clusters at different temperatures is of importance. But perhaps the main interest is to gain more fundamental understanding of the melting phenomenon, not only for finite systems, but also for bulk. Fundamental theories about cluster melting were developed in, e.g., refs. 1–6 and a large number of simulation studies have confirmed the theories.^7–13 One of the striking differences between cluster and bulk is that clusters exhibit a melting temperature interval obeying the rule that the smaller the cluster, the longer the interval.¹⁴ In this interval, the solid and molten phases coexist. The coexistence might be described by an equilibrium constant of the ratio between the number of clusters being solid and molten, respectively, in an ensemble. Small clusters are entirely solid or molten at a given moment while larger ones, as for instance Au₄₅₉, approach the behavior of bulk, i.e., a fraction of a single cluster might be molten and a fraction might be solid simultaneously.¹⁵ As the most mobile atoms in the cluster are the ones at the surface, the melting usually starts at the surface. Whether surface melting is a phase transition separate from the overall melting has been a topic of debate (see e.g.refs. 11 and 13).

Finite melting temperature intervals have been observed in experiments as well. The experimental work that has drawn most attention is the work of the Haberland group at Universität Freiburg.^16–18 In their approach, size-selected clusters assume a desired temperature by passing a heat bath of rare gas. Subsequently, a photofragmentation pattern is recorded. The temperature of the clusters can be increased either by injecting a known amount of photon energy after the heat bath or by increasing the heat bath temperature. When the photofragmentation patterns match, the photon energy corresponds to the temperature increase and the heat capacity can be calculated. In previous publications,^19–21 we verified by simulation the assumption that the clusters collide sufficiently many times with the rare gas atoms to adopt the heat bath temperature. The caloric curves of sodium clusters that were obtained by the Haberland group confirm that solid and molten cluster phases coexist within a finite temperature interval. Accurate simulations of sodium clusters are however difficult to perform since a good model potential is still missing.

Interestingly, many different methods have been used for simulation of clusters, like molecular dynamics (MD) with fixed energy,¹⁴ MD with Nosé–Hoover thermostat or a rare gas heat bath,²² Monte Carlo (MC) with fixed temperature²³ and fixed total energy,^24,25 j-walking MC⁸ and parallel tempering MC.⁹ Also, molecular dynamics simulations with the internuclear forces calculated using density functional theory have now been reported.²⁶ Furthermore, simulations aiming to obtain the free energy or the density of states like the finite-time variation method¹⁰ are also found in the literature. A majority of the publications are on Lennard-Jones clusters. The most frequently investigated cluster sizes are the icosahedrally shaped ones, but also other, e.g. cuboctahedral, clusters have been simulated.²⁷

The most common way to simulate thermal properties of an ensemble of clusters is to simulate one cluster which is isolated and prevented from fragmentation (evaporation). However, Rytkönen et al.^11,28 have shown that in a true picture of the thermal behavior of clusters, it is necessary to take into account the adsorption and desorption of cluster atoms, not only in principle, but also in practice for some clusters. Rytkönen et al. found in simulations that argon clusters with more than 309 atoms start to fragment already below the estimated melting temperature. As a solution to the problem, they performed simulations of clusters in equilibrium with its vapor. They did however find that for materials with a very low vapor pressure, such more complicated simulations are not necessary. Support for simulations of metal clusters without vapor is found in the recent publication by Moseler and Nordiek⁸ They simulated clusters with the pairwise atomic interaction represented by the Morse potential and found that the shorter the range of the potential the smaller was the temperature difference between melting and boiling. The noble gases have short-range interactions and the trend can be exemplified by the fact that the melting and boiling temperatures of bulk xenon are very close to each other (161 K and 165 K, respectively). Metal clusters, on the other hand, having a relatively large interaction range, are stable in the molten phase without significant desorption of atoms in a rather broad temperature interval.

In this work we have used the MC method to simulate a palladium cluster with a fixed temperature (canonical) or a fixed total energy (microcanonical). In the simulations, we have calculated properties such as the potential energy at various temperatures and total energies, respectively. In both ensembles, the potential energy as well as other properties, exhibit significant changes upon melting and the temperature interval and the total energy interval, respectively, where melting occurs can be identified.

The fragmentation problem has been treated differently in the two ensembles. In the microcanonical ensemble, fragmentation is impossible at total energies lower than the removal energy of one atom. This energy is however lower than the energy required for melting. Thus, fragmentation of a solid cluster is in principle possible. In our simulations fragmentation of the clusters has only been observed at energies well above the melting energy, although the simulations of the molten clusters have lasted sufficiently long to accurately probe the phase space of the molten phase. In the calculations and the results of this paper, we have only included simulations at energies lower than that of the lowest-energy simulation that led to fragmentation. In the canonical simulation fragmentation is an inevitable result at all temperatures above zero if the simulation proceeds sufficiently long. One measure to prevent fragmentation is to confine the cluster in a sphere. Such a method may however favor a spherical shape of the cluster. Here we present a method to prevent the cluster from fragmentation without favoring spherical shape. Only cluster configurations where it is possible to move from any atom to any other atom via a chain of connected atoms are regarded as non-fragmented. Any fragmented configuration is rejected. This method works well but is perhaps unnecessary for the strongly bound palladium clusters. Fragmentation attempts occurred only at temperatures well above the melting temperature.

One aim of this work is to study the similarities and differences between canonical and microcanonical MC techniques. For an infinite system the two ensembles are identical in the sence that the temperature corresponds to an exact total energy. For a finite system in the canonical system, however, the total energy is distributed and the results from the different ensembles cannot be compared directly. But using the density of states (DOS), microcanonical data may be converted into temperature-dependent data. The DOS is not known for the systems that we have simulated, but we have used a two-state model for the DOS that accurately reproduces the caloric curves of the palladium clusters. When using MC simulation of small clusters, both Gallego et al.²⁴ and Grimson²⁵ found differences between the canonical and microcanonical ensembles in the melting region. In our most accurate simulations (Pd₁₃) we have not found any evidance for differences between the ensembles.

Our simulation results are also compared to those obtained in a previous work by our colleagues Grönbeck et al.²² They investigated the melting point of palladium clusters using MD simulations with the Nosé–Hoover thermostat.²⁹ In relation to their earlier work we shall not only be able to introduce interesting methodological alternatives but also extend both the set of clusters studied and the precision of the simulations so that new results can be obtained. In this way we shall be able to partially support and otherwise modify the picture of the properties of small palladium clusters in the vicinity of the melting/freezing transition obtained in their early exploratory study.

Although simulations are successful in describing the qualitative and quantitative thermal behavior of the clusters outside the coexistence region, it has become evident that the free energy barrier between the phases may prevent an accurate determination of the melting temperature.^30,31 A second focus of this work is the character of the free energy barrier between the solid and molten phases. This barrier is a problem of considerable importance, not only for simulations of melting in clusters, but also for simulations of melting in bulk^32,33 and in general for simulations of any systems with several comparable deep local free energy minima.

The article is organized as follows: In Sec. 2 we present the computational and simulation details. In Sec. 3 we discuss a model for the density of states which we use for qualitative and quantitative interpretation of the simulation results. The conversion of data from microcanonical simulation to temperature-dependent data is derived in Sec. 4. In Sec. 5 we present the results for the icosahedral clusters with a discussion of the free energy phase barriers. Results for the non-symmetrical clusters are given in Sec. 6 and conclusions are drawn in Sec. 7.

II. Monte Carlo simulation of clusters

The system of our concern is a low-pressure gas of single-sized palladium clusters. Such a gas would not be stable in a chamber since when clusters collide with each other or with the chamber wall, they readily exchange atoms so that the monodisperse nature of the gas is lost. In contrast, free clusters are experimentally mostly studied in cluster beams in vacuum. In the supersonic beam often used in cluster experiments,³⁴ the translational and rotational temperatures are considerably cooled.^35,36 Thus the vibrations are not in equilibrium with the rest of the degrees of freedom. Still, the cluster gas is a convenient system to simulate and since the translational motion is independent of that of the internal degrees of freedom, simulation results may be converted from one to other translational temperatures. Furthermore, one may expect the rotational energy to couple weakly with the vibrations. The simulations are concerned with the vibrations of the cluster. Thus the simulation results at one vibrational temperature could be combined with any translational temperature and rotational temperature. The results presented below refer to the cluster, inclusive of translation and rotation, with totally equilibrated degrees of freedom.

Since the pressure is assumed to be very low, the clusters encounter each other very rarely. Each cluster can therefore be regarded as a microcanonical system. As an alternative to the study of a canonical ensemble of clusters, it might be of equal relevance to study the properties of a microcanonical ensemble of single clusters. However, the temperature dependence of properties is usually more relevant than energy dependence in experiments. Hence, a conversion of microcanonical simulation data to canonical form is appropriate for practical use. In this paper the conversion is performed to test the correspondence of the two ensembles.

A. Monte Carlo simulation of clusters at fixed temperature

By means of simulation, we have investigated how properties, as for instance the potential energy, of free clusters with N atoms per cluster evolve with the temperature T where T varies from 50 K to 2000 K. Since the clusters are at a very low pressure, we can assume that they never collide with each other. Thus, each cluster is an isolated system with a fixed total energy. In the canonical ensemble of such clusters, the energy is thermally distributed. Let us now use one cluster to represent the whole ensemble of clusters. Then the spatial configurations of the cluster, R = (x₁,y₁,z₁,x₂,…,z_N), represent the distribution for the whole canonical ensemble. The distribution is given by the Boltzmann factor,

probability(R) ∝ exp(−U(R)/k_BT), (1)

where U is the potential energy of the single cluster and k_B is Boltzmann's constant. The momenta and kinetic energy need not be considered as the potential and kinetic energies are separable in the canonical ensemble and the momenta will contribute to the energy simply by E_kin = 3Nk_BT/2.

Besides the temperature and the number of atoms, the volume of the system should be fixed in the canonical ensemble. As the cluster gas is assumed to be ideal, the volume per cluster should be large or even infinite. A fundamental problem is that in an infinite volume, the equilibrium state of the single cluster described by eqn. 1 corresponds to all cluster atoms having evaporated from the cluster for all T > 0 (cf.refs. 11 and 28). The reason is that although the Boltzmann factor might be considerably larger for the minimum energy configuration than for fragmented configurations, the volume of the phase space corresponding to fragmented configurations is infinite. In order to avoid the inevitable fragmentation, the standard Metropolis Monte Carlo method³⁷ therefore has to be modified so that a move of an atom that leads to fragmentation is rejected. Hence, we define the average value of, for instance, the potential energy at the temperature T as

where the sampling is over all configurations Rthat describe a cluster.

One method to prevent one atom from leaving the cluster is to reject a configuration if any atom is more distant from its nearest neighbor than a user-defined maximal distance. This method does however not prevent an escape of e.g. a dimer from the rest of the cluster. One method to define a cluster that prevents escape of any atom is to require that all atoms have to be within a sphere with fixed radius from the centre of mass.^38,39 But with such a definition of a cluster, the shape of the cluster is somewhat constrained, favoring spherical shape. As will be shown in Sec. 5, the shape of molten clusters tend to be non-spherical. Therefore, we have used an alternative definition of a non-fragmented cluster by which the cluster is constrained to be connected but no constraint is put on cluster shape: A configuration of atoms form a cluster if all atoms are linked to all other atoms via a chain of adjacent atoms. Two atoms which are closer to each other than an element-specific limit distance are said to be adjacent. To determine whether a configuration is connected is exactly the problem described in graph theory (see e.g.ref. 40). Our algorithm for separating fragmented configurations from non-fragmented ones is described in the appendix.

The disadvantage of the canonical ensemble is that we have to introduce some kind of extra parameter to distinguish between a cluster and a fragmented cluster. In our case the parameter is the interatomic limit distance for two atoms to be adjacent. However, the interaction between the palladium atoms is so strong that the process to fragment a cluster is very slow due to the low vapor pressure of the cluster.¹¹ Hence, we have seen fragmentation attempts on very few occasions only. In our study, Pd₁₃ is the cluster for which an attempt to fragment occurred at the lowest temperature, which was at 1500 K. This temperature is well above the melting temperature of the cluster. In the appendix we find that a limit of 4.0 Å is a good divisor for fragmented and non-fragmented clusters.

B. Monte Carlo simulation of a cluster at fixed total energy

In the microcanonical MC simulation, the system to simulate is a cluster which is isolated in an infinite volume. The total energy of the cluster is fixed at E. For the sampling of quantities we have used the efficient microcanonical sampling (EMS) method developed by Severin et al.⁴¹ A detailed description of this method is also found in ref. 42. A short derivation of the method is given here.

One of the fundamental principles of statistical mechanics is that for a system at fixed total energy, all states with this total energy are equally probable. A state is here defined by the spatial coordinates of the cluster atoms R and the momenta of the atoms P = (p_x1,p_y1,p_z1,p_x2,…,p_zN). However, no quantities related to the momenta, except the kinetic energy, are of interest in this study. Since the potential energy is sampled and the total energy is fixed, the average kinetic energy is automatically obtained and there is no need for simulating the momenta. At a total energy E, the probability of a spatial configuration R, regardless of momenta, must be proportional to the number of available momentum configurations P, on the condition that the kinetic energy satisfies E_K = E − U. The number of momentum configurations is proportional to the density of momentum states,

Ω_K(E_K) = constant × E^3N/2−1_K. (3)

Thus, the probability of a spatial configuration R satisfies

probability(R) ∝ (E − U(R))^3N/2−1. (4)

Standard Metropolis Monte Carlo simulation³⁷ can now be used. One atom at a time is moved a small distance and the new spatial configuration is accepted if where S[0,1] is a uniform random number between zero and one.

C. The many-body alloy potential

The interaction between the atoms in the cluster is described by the many-body alloy potential, developed by Tománek and co-workers.⁴³ The potential energy of the cluster is written as a sum over all atoms in the clusterwhereHere r_ij is the distance between atoms i and j. The repulsive part of the binding energy is considered to be pairwise and of Born–Mayer type. The attractive part is modeled by a tight-binding second-moment approximation of the d band density of states. This model leads to a many-body character of the attractive part of the interaction. The parameters that we use were calculated by Tománek et al.⁴³ by fitting the energy to LDA calculations of totally frozen bulk palladium. The parameter values are ξ₀ = 1.2630 eV, ε₀ = 0.08376 eV, r₀ = 2.758 Å, p = 14.8 and q = 3.40. These parameter values were also used in the work by Grönbeck et al.²² with which we compare our results. This potential compares favorably with SCF-CI calculations of the optimized geometry⁴⁴ for finite systems like Pd₁₃. The mean radius of the icosahedral structure of Pd₁₃ is with SCF-CI 2.60 Å and with the many-body alloy potential 2.50 Å. Moreover, the energy difference between the icosahedral and octahedral isomers is 0.12 eV with both techniques.

It must be mentioned that other sets of parameters for this potential and palladium have been used by others. For instance, Rey and coworkers⁴⁵ used other parameters and therefore obtained another melting temperature (816 K) for Pd₁₃ to be compared with 1110 K for Pd₁₃ with our parameters.

D. Simulation details

1. Simulation at fixed temperature. Three kinds of simulation have been performed. The first kind was carried out at fixed temperature and is referred to as MC-CAN. They are run in the following way: the cluster is prepared in a configuration which is equal or close to the minimum-energy configuration. For the icosahedral clusters, Pd₁₃, Pd₅₅ and Pd₁₄₇, this is accomplished by starting in a near-icosahedral structure and subsequently finding the minimum in potential energy by using a steepest-descent method. For Pd₃₄ and Pd₅₄, there are probably many configurations with an energy close to the global minimum. By an annealing procedure, the molten clusters were brought to 1 K and then a steepest-descent search was started. For Pd₅₄, the minimum-energy configuration that was found is simply the icosahedral structure of Pd₅₅ with one outer atom removed. For Pd₃₄, the minimum-energy configuration that we found has five core atoms and one surrounding non-spherical shell with 29 atoms.

A simulation is started by letting the minimum-energy configuration approach equilibrium at desired temperature. New configurations are generated by letting the atoms move randomly one by one with a small maximal displacement. During the meltdown stage of the simulation when no statistics are collected, an algorithm is used to adjust the maximal displacement aiming at finding a value which will lead to a rejection of approximately half the generated configurations during the succeeding sampling according to standard Monte Carlo technique.⁴⁶ The maximal displacement was changed ten times during the meltdown and thereafter held constant in the main part of the simulation. In addition to the rejection due to a small Boltzmann factor, trial configurations are rejected if they form fragmented clusters in accordance to Sec. 2.A. The meltdown part of the simulation was run for about 1 million to 3 million single-atom move attempts (“steps”) depending on cluster size.

After the meltdown, quantities such as the potential energy are sampled for a number of steps. During the sampling period the maximal displacement was kept at a fixed value based on the rejection rates obtained during meltdown. The resulting rejection rate during the sampling period was confirmed to fall between 48% and 52% in all the simulations. The number of steps during the sampling varied normally from 10 million to 30 million. In addition to the potential energy, U, the mean distance to the centre of mass, the radii of gyration, the specific heat C_v and the Lindemann index δ were sampled. The mean distance to the centre of mass is defined as

where x_i is the position of atom i and the centre of mass is x_CM. This quantity can describe a volume expansion of the cluster. The radii of gyration are calculated for the principal axes and are defined as where I_i is the moment of inertia along principal axis i. Note that if all the atoms in the cluster have the same mass, then the total mass of the cluster, m, is included also in I_i and is canceled in r_i. For a solid ellipsoid, the relation between the semi-axes a, b and c and the radii of gyration is Thus, if approximating the cluster with an ellipsoid, the radii of gyration can give information about the shape of the cluster. The quantity δ measures the fluctuations in distances between atoms and is defined as If any atoms switch positions in the cluster, δ will increase drastically and thus a rise in δ is an indicator of melting.

2. Simulation at fixed total energy. The second kind of simulation was carried out at fixed total energy and is referred to as MC-MICRO. They are run in analogy with the MC-CAN simulations described above. The initial configuration of the clusters were the same as for MC-CAN, i.e., the configurations corresponded approximately to the minimum energy state. In the MC-MICRO simulations the icosahedral cluster Pd₃₀₉ was included. A meltdown preceded the sampling of properties. The meltdown part was run for about 1 million to 15 million steps depending on cluster size. The number of steps during the sampling varied from 10 million to 50 million. In the coexistence region, some simulations lasted even longer. In the simulations the potential energy, the mean distance to the centre of mass, the radii of gyration and the Lindemann index were sampled. In analogy with the MC-CAN simulations the maximal displacement was adjusted during the meltdown and a fixed maximal displacement was used during the sampling giving approximately 50% rejection of generated configurations.

3. Simulation with annealing temperature. The third kind of simulation was carried out at a constantly decreasing temperature and is referred to as MC-ANNEAL. These simulations are similar to those denoted MC-CAN. However, during the sampling period the temperature is decreased before each MC step by ΔT = (T_final − T_initial)/N_s, where N_s is the number of steps during the sampling. Moreover, the maximal displacement is decreased nine times during the sampling period in order to avoid a too high rejection rate when the temperature is lowered. If ΔT is large, the cluster will transform from a molten phase to a glass-like phase. If ΔT is small enough the cluster might attain the highly symmetric solid structure. The phase switch is however obstructed or even prevented by a free energy barrier.

III. A two-state model for the density of states

In an infinite system, the transition from one phase to another at a given pressure, occurs at a specific temperature, barring any metastabilities. In contrast, for a finite system like the cluster, the solid (s) and molten (m) phases coexist within a certain temperature interval. Jellinek et al.¹⁴ concluded that in simulations of Ar₁₃, the time spent between every switch of phase was long enough in order to equilibrate in that phase. Hence we can speak of switches between fully developed phases.

The equilibrium between the phases at a fixed temperature can be thought of as a unimolecular reaction equilibrium

Pd_N(s) ⇌ Pd_N(m), (14)

where the equilibrium constant isHere F_m = U_m − TS_m is the Helmholtz free energy of one molten cluster. The partition functions are for single clusters and their form in classical mechanics isHere Ω_s and Ω_m are the densities of states (DOS) of the solid and molten isomers, respectively. E₀ is the minimal attainable (potential) energy of the cluster. Although the melting occurs over a temperature interval, we define the melting temperature, T_M, as the temperature where the equilibrium constant satisfies K_T = 1. When the number of atoms increases, the DOS becomes enormous and K_T switches from 0 to ∞ in an infinitesimal temperature interval, in accord with thermodynamics.

An equilibrium between coexisting phases may be established at a fixed total energy as well. The microcanonical equilibrium constant is given by

At a certain total energy the DOS of the two phases are equal. This energy we call the melting energy, E_M.

In previous publications^19–21 based on the caloric curves of Grönbeck et al.,²² we proposed the following model for the DOS of the icosahedral clusters Pd₁₃ and Pd₅₅:

Ω(E) = Ω_s(E) + Ω_m(E) ∝ (E − E₀)^n₁ + k₂(E − E₀)^n₂, (19)

where E is the total energy of the cluster. The model is motivated by the assumption that the cluster might be in either of two states, one entirely solid structure or one of many entirely molten structures which are assumed to be statistically equivalent. The DOS is then expressed as the sum of independent contributions for the two states. The vibration, translation and rotation of the cluster give rise to power law terms like (E − E₀)ⁿ, but to account for the anharmonicity of the vibrations, the exponents n₁ and n₂ are adjustable parameters. Although E₀ may be calculated by a search for the ground state configuration of a cluster, it is regarded an adjustable parameter as well. The prefactor k₂ is the last parameter (the prefactor of the solid phase term is arbitrarily set to unity).

A. The temperature-dependent caloric curve

The temperature-dependent caloric equation based on this model will be where Γ(x) is the gamma function. At low temperatures, 〈E〉_T − E₀ ∼ (n₁ + 1)k_BT, and for high temperatures, 〈E〉_T − E₀ ∼ (n₂ + 1)k_BT. In-between, the caloric curve of the model DOS exhibits the typical S-shape of melting as seen in the upper panel in Fig. 1. At low temperatures the vibrations in the cluster are harmonic and 〈E〉_T − E₀ ∼ (3N − 3)k_BT. However, with increasing T, the atoms start to probe the anharmonicity of the potential and n₁ + 1 will generally be greater than 3N − 3. This is even more true for the molten phase and n₂/n₁ is about 1.3 for Pd₁₃, Pd₅₄, Pd₅₅, Pd₁₄₇ and Pd₃₀₉. The fitting parameters in the model DOS are obtained by a least-squares fit of simulation data to the expression in eqn. (20).

Fig. 1 Upper panel: The total energy per atom for Pd₁₃ (rings), Pd₅₅ (stars) and Pd₁₄₇ (triangles) versusT according to the MC-CAN simulations. The plusses represent the outcome from the MD-NH simulations in ref. 22. The solid lines are the caloric curves given by eqn. (20). The parameters of the model DOS are adjusted to fit the MC-CAN energy data. The vertical lines indicate the melting points given by eqn. (21). From the left, the lines correspond to Pd₅₅, Pd₁₄₇ and Pd₁₃, respectively. Lower panel: C_v per atom from the MC-CAN simulations of Pd₁₃ (rings), Pd₅₅ (stars) and Pd₁₄₇ (triangles) versusT. The solid lines are derived from eqn. (20). NB: logarithmic scale on the ordinate.

Using the model DOS in eqn. (19), the melting temperature giving K_T = 1 is

An alternative way of defining the melting temperature is to use the temperature for which the heat capacity is maximal. The heat capacity is the derivative of eqn. (20),For a finite system like a cluster, C_v does not show a singularity at the melting temperature, but rather a peak of finite width. Employing the model DOS (eqn. (19)), the maximal heat capacity is obtained when The difference in melting point determined by eqns. (21) and (23), respectively, due to the last factor is found to be negligible. One of the advantages using the canonical ensemble is that the heat capacity can be obtained directly from a simulation since it is proportional to the fluctuation in energy:Eqn. (24) is used for calculating the heat capacities in the lower panel of Fig. 1 as it is more accurate than using numerical differentiation of the energy in the upper panel of Fig. 1.

B. The energy-dependent caloric curve

The caloric curve in the microcanonical ensemble is the kinetic energy versus total energy. In order to optimize the parameters we need an analytical expression for 〈E_K〉_E in the terms of the model. The kinetic energy is where f(E_K|E) is the distribution of E_K on the condition that the total energy is E. E − E₀ is the available energy that should be distributed among kinetic and potential energies. f(E_K|E) may be expressed as an integral over the total coordinate phase space and the total momentum phase space:In this section δ is Dirac's delta function. The denominator is proportional to the DOS of the total energy as,Here h is Planck's constant. Let us define the density of states of the kinetic energy and the potential energy⁴⁷ in analogy with eqn. (27),The constant h_K should have the dimension (mass length time⁻¹) and h_p the dimension (length). In this derivation the constants h_K and h_p will always appear as h_Kh_p and the values of the constants may be arbitrary chosen as long as h_Kh_p = h. With these definitions, the distribution function in eqn. (26) can be written asFurthermore, these definitions satisfy the convolutionAccording to the Laplace transform of a convolution we have

[script L][Ω_kin(E_K)][script L][Ω_pot(U)] = [script L][Ω(E)] (32)

Since there are 3N momenta contributing to the kinetic energy, the DOS of the kinetic energy is

Ω_kin (E_k) ∝ E_k^3N/2−1 (33)

Using the model DOS for Ω(E) in eqn. (19) and the Laplace transform, the model DOS of the potential energy is obtained as Now we have an analytical expression for f(E_K|E) in terms of the model and finally, the caloric curve in the microcanonical ensemble according to the model DOS is found to beThis model gives a graph of the kinetic energy versus the total energy exhibiting three parts. At low and high total energy, the cluster is entirely solid and entirely molten, respectively. The curve is then linear according to For intermediate energies the two phases of the cluster coexist and the curve performs a wiggle. The expression in eqn. (35) can be fitted to simulation data determining the parameters of the model. When the parameters are calculated, the melting energy according to the model DOS can be calculated from

K_m(E_M) = 1 ⇔ Ω_m(E_M) = Ω_s(E_M) ⇔ (E_M − E₀)^n₁ = k₂(E_M − E₀)^n₂ ⇔ E_M = E₀ + (k₂)^{1/(n₁−n₂)} (37)

IV. Conversion of energy-dependent data into temperature-dependent data

For large systems, the distribution of kinetic energy at constant total energy, as well as in the canonical ensemble, is extremely narrow. The canonical and microcanonical ensembles are then equivalent and the total energy can be immediately translated into temperature. In contrast, in the case of a smaller system as Pd₁₃, the fluctuation in kinetic energy is considerable and a Boltzmann weighting must be employed to convert microcanonical data into canonical data. For such a Boltzmann weighting, the DOS of the cluster must be known. In the microcanonical simulations presented here, no direct determination of the DOS is included but by using the fitted two-state DOS model, the conversion of microcanonical data into canonical data can be carried out.

In the microcanonical simulations, quantities are averaged at a constant total energy. Say that A is a quantity that has been sampled for energies up to E = E_max. The average of A at the temperature T is then given by a Boltzmann weighting including the density of states, i.e.

The approximation of replacing ∞ with E_max is acceptable if the temperature is so low that the contribution from E > E_max is negligible. We have chosen to define the highest temperature T_max for which the conversion is acceptable as the temperature givingUsing data from the microcanonical simulations and the model DOS with the fitted parameters, T_max can be calculated and subsequently 〈A〉_T for T < T_max.

V. Results for Pd₁₃, Pd₅₅, Pd₁₄₇ and Pd₃₀₉

Let us first have a look at the caloric curve from MC-CAN simulations of Pd₁₃, Pd₅₅ and Pd₁₄₇. The caloric curve of the total energy is constructed by adding E_kin = 3Nk_BT/2 to the sampled potential energy and it is drawn as big symbols in the upper panel of Fig. 1. As expected, the curves exhibit the typical S-bend around the melting point. Simulation results are fitted to eqn. (20) in a least-square sense and the resulting parameters are given in Table 1. The fitted model caloric curves are also shown as solid lines in Fig. 1. The curves are in good agreement with simulation data and we find that the model DOS describes the thermal behavior of the clusters very well. As a measure of the accuracy of the DOS model, the maximal difference between the potential energy from the simulations and from the model DOS,

Q_Etot = max|E_{tot,simulation} − E_{tot,model DOS}|, (40)

was calculated and given in Table 1. The melting points in accordance with eqn. (21) are indicated by vertical lines.

Table 1 The model DOS parameters and the melting points after fitting MC-CAN simulation data to eqn. (20). (Exception: Pd₃₄.) As a comparison, the minimum energy per atom obtained by a steepest descent search is listed as E_min. The melting points of the MD-NH simulations²² are obtained using eqn. 41

N	n1/N	n2/N	k2/eV^(n1−n2)	E0/N/eV	(QEtot/N)/eV	Tm/K	Tm,MD-NH/K	Emin/eV
13	3.10	4.14	2 × 10^−9	−2.67	0.011	1110	1099	−2.66
55	3.46	4.47	7 × 10^−69	−2.99	0.022	908	1034	−2.97
147	3.42	4.57	3 × 10^−293	−3.12	0.016	1076	1163	−3.11
34	—	—	—	—		≈350	—	−2.85
54	3.24	4.33	3 × 10^−71	−2.97	0.013	893	—	−2.95

---

The results from ref. 22, where Grönbeck et al. employed molecular dynamics simulation with the Nosé–Hoover thermostat (MD-NH), are drawn as plusses. The data from ref. 22 is slightly modified as compared to Fig. 1 in ref. 22. In that publication, the temperature of the cluster was calculated from E_kin = 3Nk_BT/2 despite the fact that the rotation and translation were kept to zero. Here, we use the correct expression

E_kin = (3N − 6)k_BT/2, (41)

when setting the temperature scale of the MD-NH simulations and by doing so we find good agreement between the MC-CAN and the MD-NH simulations. Precisely at melting, the MC and MD-NH curves do however disagree for the two larger clusters. This is not due to fundamental differences between the methods, but is a result of the simulations in the MD-NH case being too short. The temperatures for which the caloric curves differ are within the coexistence region. There is a free energy barrier between solid and molten structures and if the cluster is large the simulations might need to be unreasonably long to allow for a large number of phase switches and a correct distribution between the phases. As explained in detail below we argue that the barrier does not prevent a correct simulation in the coexistence region in the case of Pd₁₃ but does so in our simulations of Pd₅₅ and Pd₁₄₇. Since the simulations start with the solid structure, this structure was favored leading to overestimated melting points. Consequently, the longer MC-CAN simulations lead to lower melting points for Pd₅₅ and Pd₁₄₇. The melting points from the MD-NH simulations obtained using eqn. (41) are listed in Table 1. For Pd₁₃ there is good agreement with the MC-CAN simulations. For Pd₅₅ and Pd₁₄₇ the melting points are probably overestimated in both the MC-CAN and MD-NH simulations.

The heat capacity per atom is obtained as the variance of the energy during a simulation (see eqn. (24)). In the lower panel of Fig. 1, C_v per atom from the MC simulations is compared with the heat capacity derived from the model DOS. (Note the logarithmic scale of the heat capacity). The agreement in determination of the melting point is good. The peaks in heat capacity are due to the fluctuations between the solid and molten phases that result in fluctuations in potential energy. The heat capacity from simulation and from the DOS model differ qualitatively at low temperatures. The model DOS implies that the heat capacity should be constant at low temperatures. In the simulations, however, the heat capacity increases with temperature as the vibrations in the cluster become increasingly anharmonic when the temperature is raised.

The caloric curves from the microcanonical simulations MC-MICRO are drawn in Fig. 2. The upper curve in each panel is the potential energy minus the minimal potential energy and the lower curve is the kinetic energy of the clusters. As predicted by the model DOS, the curves show two regions of linear character and an intermediate coexistence region. At the melting, the kinetic energy behaves peculiarly in the sense that it decreases when the total energy increases. This is a well-known phenomenon for finite systems.^3,4 The explanation is that in principle only two states are available. At the switch from solid to molten phase, energy must be taken from the kinetic energy and given to the potential energy in order to attain the high-energy molten structure. The kinetic energy data in Fig. 2 are fitted to eqn. (35) giving the parameters in Table 2. The agreement with the parameters from MC-CAN and MC-MICRO is striking and so is the calculated melting temperature, except for Pd₁₄₇. The caloric curves calculated using eqn. (35) are drawn as solid lines in Fig. 2. The agreement is good showing that the two-state model DOS indeed captures the fundamental behavior of the clusters. In Sec. 5.A, we will return to the discrepancy for the parameters of Pd₁₄₇ which is due to the barrier between the solid and molten phases and not to the character of the ensembles. As a measure of the accuracy of the DOS model, the maximal difference between the potential energy from the simulations and from the model DOS,

Q_Ekin = max|E_{kin,simulation} − E_{kin,model DOS}|, (42)

has been calculated and is given Table 2.

Fig. 2 The caloric curves of the microcanonical simulations for Pd₁₃, Pd₅₅, Pd₁₄₇ and Pd₃₀₉. The rings show the kinetic energy and the pentagrams show the potential energy minus E₀versus total energy. The solid lines are the caloric curves resulting from the model DOS according to eqn. (35). In the figure of Pd₃₀₉, grey lines show the caloric curves assuming that the cluster melts bulk-like. The inset is a magnification of the melting interval of Pd₃₀₉.

Table 2 The model DOS parameters and the melting points after fitting MC-MICRO simulation data to eqn. (35). (Exception: Pd₃₄.) As a comparison, the minimum energy per atom obtained by a steepest descent search is listed as E_min. The melting point T_M is calculated as in eqn. (21)

N	n1/N	n2/N	k2/eV^(n1−n2)	(E0/N)/eV	(QEkin/N)/eV	TM/K	(Em/N)/eV	[(Em − E0)/N]/eV	Minimum energy per atom/eV
13	3.14	4.14	3 × 10^−9	−2.67	0.005	1116	−2.32	0.35	−2.66
34	—	—	—	—	—	—	−2.76	0.09	−2.85
54	3.08	4.07	3 × 10^−63	−2.96	0.007	894	−2.68	0.28	−2.95
55	3.28	4.32	5 × 10^−70	−2.98	0.006	926	−2.67	0.30	−2.97
147	3.42	4.59	7 × 10^−293	−3.12	0.015	996	−2.78	0.34	−3.11
309	3.31	4.45	3 × 10^−727	−3.17	0.020	1121	−2.79	0.37	−3.17

---

The correspondence between the canonical and microcanonical ensembles can be seen in Fig. 3. The potential energy in the MC-CAN simulations are shown as symbols and the converted data from MC-MICRO simulations according to eqn. (38) are shown as solid lines. The agreement is good for Pd₁₃ and Pd₅₅. Again, the discrepancy for Pd₁₄₇ is due to the phase barrier.

Fig. 3 The simulated average potential energy is shown as a function of temperature. The solid lines show the converted microcanonical data (eqn. (38)). The lines are from the bottom: Pd₃₀₉, Pd₁₄₇, Pd₅₅ and Pd₁₃. Simulation data from the canonical MC-CAN runs are shown as rings (Pd₁₃), pentagrams (Pd₅₅) and squares (Pd₁₄₇).

The geometric shape of a cluster changes upon melting.¹⁴ In Fig. 4, the average radii of gyration versus temperature are drawn for Pd₁₃, Pd₅₅, Pd₁₄₇ and Pd₃₀₉. The radii of gyration are calculated for each configuration and are then sorted in order to sample the maximal, medium and minimal radius of gyration. In the solid clusters, the radii of gyration are almost identical. In contrast, the radii of gyration for molten clusters are split up with one radius even decreasing compared to the solid phase and the other ones dramatically increasing. The radii of gyration of Pd₁₃ at 1300 K are in ratios 2.5 ∶ 2.3 ∶ 2.0, and if approximating the cluster with a solid ellipsoid, the corresponding semi-axes are 4.4 ∶ 3.5 ∶ 2.8, which is far from a spherical shape. The semi-axes of Pd₁₄₇ at the same temperature are in ratios 9.5 ∶ 8.3 ∶ 7.5. The transformation from spherically symmetric shapes to non-spherical shapes upon melting contradicts the general knowledge that a liquid droplet favors a spherical shape due to surface tension.⁴⁸ However, for these tiny clusters, the configuration fluctuations are noticeable, causing a split up of the radii of gyration. We believe that the ratio of the fluctuations to the cluster radius will diminish when the number of atoms increases so that the cluster eventually behaves like a liquid droplet. In Fig. 4, the ratio of the maximal to the minimal radius is indeed smaller for Pd₁₄₇ than for Pd₁₃. The rings in the Pd₅₅ figure are from the MD-NH simulations in.²² The two methods agree except in the temperature interval 900–1050 K, where the two methods also disagree in energy (Fig. 1). The reason for the disagreement is also the same, i.e., the cluster spends too short a time in the molten phase in the MD-NH simulations which are not sufficiently long.

Fig. 4 The maximal, medium and minimal radii of gyration versusT for Pd₁₃, Pd₅₅, Pd₁₄₇ and Pd₃₀₉. The star-lines show canonical MC-CAN results. The thick grey lines are the converted microcanonical MC-MICRO results. The vertical lines indicate the melting points in Table 1 (for Pd₃₀₉ according to Table 2). The curves with rings for Pd₅₅ are from the MD-NH simulations in ref. 22.

In Figs. 5 and 6, δ, the dispersion in atomic separation, and d_CM, the mean distance to the centre of mass, are drawn. It is noticeable that the rise in δ for Pd₁₃ starts far below the melting point. The rise in δ starts when surface atoms start to switch positions and this occurs at the very beginning of the melting interval (800 K). The dramatic change in shape does however not occur until at 1000 K. The dramatic change in shape of Pd₁₄₇ upon melting is depicted in Fig. 7. Two snapshots are shown for 100 K (solid phase) and 1400 K (molten phase), respectively.

Fig. 5 The Lindemann index, δ, versusT for Pd₁₃, Pd₅₅, Pd₁₄₇ and Pd₃₀₉. The star-lines show canonical MC-CAN results. The thick grey lines are the converted microcanonical MC-MICRO results. The dashed grey line for Pd₁₃ shows δ when calculated from converted microcanonical 〈r_ij〉 and 〈r²_ij〉. The vertical lines indicate the melting points in Table 1 (for Pd₃₀₉ according to Table 2).

Fig. 6 The average distance to the centre of mass, d_CM, versusT for Pd₁₃, Pd₅₅, Pd₁₄₇ and Pd₃₀₉. The star-lines show canonical MC-CAN results. The thick grey lines are the converted microcanonical MC-MICRO results. The vertical lines indicate the melting points in Table 1 (for Pd₃₀₉ according to Table 2).

Fig. 7 Two snapshots of Pd₁₄₇. The right one is at 100 K (solid phase) and the left one at 1400 K (molten phase).

For the radii of gyration, δ and d_CM, the agreement between the MC-CAN results and the converted MC-MICRO results is good with two exceptions. As already mentioned, the MC-CAN and MC-MICRO simulations disagree for Pd₁₄₇ in the melting region due to the phase barrier. The second discrepancy for δ is found for Pd₁₃ and Pd₅₅ near melting. In both cases the rise is at lower temperatures with the canonical data. However, δ has the nature of a standard deviation. In the MC-MICRO simulations, this standard deviation is calculated for subsets of the phase space of the cluster and then averaged into the temperature-dependent values. Such a procedure cannot be expected to give a correct correspondence with the canonical value in the melting region. A better way would be to convert the averages 〈r_ij〉_E and 〈r²_ij〉_E into temperature-dependent 〈r_ij〉_T and 〈r²_ij〉_T and thereafter to calculate a temperature-dependent δ. When doing so for Pd₁₃ the dashed curve in Fig. 5 is obtained. The agreement with canonical data is now better but not perfect. A further explanation of the discrepancy between canonical and converted microcanonical results is that although the simulations were sufficiently long to give accurate averages for most quantities, even longer simulations are required to give accurate values of δ. For instance, we see that the canonical δ is not even constantly increasing with temperature as must be expected. Since the phase space of the cluster is finite we believe that if the simulations were infinite, a perfact match between canonical and converted microcanonical δ would result. Also in the results from Gallego et al.²⁴ and Grimson²⁵ we can see that the agreement between canonical and microcanonical data is much better for the potential energy than for δ. We believe that a direct conversion of δ and too short simulations are the reasons for the discrepancy for δ also in their simulations.

An indication that the cluster can be decomposed into two major isomers is the distribution of the potential energy for Pd₅₅ at 900 K as shown in Fig. 8. The potential energy density is split up into two peaks. The solid line is the distribution according to the model DOS (eqn. 19). Simulation results and the model are in fair agreement but the relative heights of the two peaks are very sensitive to temperature. The model is not capable of capturing this temperature behavior exactly. A similar histogram with two peaks is found for Pd₅₅ in the MC-MICRO simulations at a total energy corresponding to the coexistence region.

Fig. 8 Distributions of the potential energy for Pd₅₅ at 900 K according to the canonical MC-CAN simulation. Averages over 1000 steps in the MC simulation are collected in the histogram. The distribution of potential energy for a single configuration according to the model DOS is drawn as a solid line.

A. Free energy barrier

In the coexistence region the cluster constantly switches between solid and molten phases. This can be seen for instance by following how the potential energy varies during a simulation. The mean distance to the centre of mass and the potential energy of Pd₁₃ are plotted versus number of macro steps in Fig. 9. The cluster temperature is 1100 K which is in the middle of the coexistence region. (One macro step equals N MC steps where N is the number of atoms in the cluster. The macro step is a better quantity when comparing simulation lengths of clusters of different sizes.) Fig. 9 clearly shows how the cluster frequently switches phases. The ratio of macro steps in molten phase to solid phase equals the equilibrium constant. However, in the cases of Pd₅₅ and Pd₁₄₇, only a few switches are found and in the case of Pd₁₄₇, we have never observed any true switch from molten to solid phase in a MC-CAN simulation (Fig. 10). This problem is due to the free energy barrier between the regions in the phase space. As a consequence of the barrier, it is difficult to determine the equilibrium constant K_T as a function of T, and particularly the determination of the melting point is inaccurate. In contrast, the barrier for switches between the phases for Pd₁₃ is low enough that, in the course of a normal simulation, the cluster can switch phases so many times that proper statistics should have been obtained. The problem to obtain correct distribution between phases is the explanation why the MD-NH simulations of Grönbeck²² differ from the MC simulations in the coexistence region. The MD-NH simulations were interrupted too early so that the cluster spent a shorter time in the molten phase than in the MC simulations. The consequence is an estimated higher melting point. In the case of Pd₁₃ the barrier is easily overcome and the MD-NH and MC simulations agree.

Fig. 9 The mean distance to the centre of mass (upper panel) and potential energy (lower panel) of Pd₁₃ in the course of a MC-CAN simulation are shown at 1100 K. Each dot is an average over 1000 simulation steps.

Fig. 10 Upper panel: The potential energy of Pd₁₄₇ in the course of a MC-CAN simulation is shown at 1100 K. Each dot is an average over 10 000 simulation steps. Lower panel: The potential energy of Pd₅₅ in the course of a MC-CAN simulation is shown at 900 K. Each dot is an average over 1000 simulation steps.

For Pd₁₃ and Pd₅₅ the MC-CAN and MC-MICRO simulations are in good agreement in the coexistence region. This is not the case for Pd₁₄₇. Presumably the reason is the phase barrier. In Fig. 11 the evolution of the potential energy is followed in a MC-MICRO simulation at melting. Comparing the energy scales for Pd₁₄₇ in Figs. 10 and 11, we see that the phase switch implies a smaller change in potential energy in the microcanonical case than in the canonical case. At fixed total energy, the phase switch occurs between relatively high-energy solid and low-energy molten phases while at fixed temperature both low-energy solid phases as well as high-energy molten phases are attainable. The phase switch to molten phase seemingly occurs at a lower “temperature” in MC-MICRO than in MC-CAN. The reason might be the difference in the potential energy step between the phases. The reason might however also be mere chance which is supported by the fact that the melting point of Pd₅₅ is higher in the MC-MICRO simulations.

Fig. 11 The potential energy of Pd₁₄₇ in the course of a MC-MICRO simulation is shown at the total energy −2.765 eV per atom. Each dot is an average over 10 000 simulation steps.

Let us return to the caloric curve of Pd₃₀₉ in Fig. 2. The curves show the S-bend at melting which is typical for a two-state system. In bulk, the behavior is different. Since bulk solid can melt gradually so that for instance half the atoms of the bulk are in a molten phase, bulk melting in the microcanonical ensemble is characterized by a horizontal line in a graph of the kinetic energy versus the total energy. Such a line is indicated in Fig. 2 for Pd₃₀₉. The reason why Pd₃₀₉ seems to melt like a two-state system might be any of two explanations or both. Obviously, the cluster might be sufficiently small so that it is still well described as a two-state system. However, the apparent two-state behavior might be an artifact due to the phase barrier. Although the simulations of Pd₃₀₉ lasted for 50 million steps, no fluctuation between phases was ever observed. In all the simulations the starting configuration was a solid one and thus it might be so that the cluster does not manage to overcome the barrier until the total energy is high enough that the entire cluster should be molten. If this is the case and the simulations are sufficiently long so that the coexistence between fractions of solid and molten phases is established, the caloric curves of Pd₃₀₉ should be like the grey lines in Fig. 2. The conclusion is that our simulations cannot resolve whether the Pd₃₀₉ is a true two-state system or a bulk-like system. An example of a rather small cluster with a bulk-like melting is Au₄₅₉ in the simulations of Cleveland et al.¹⁵

The obstacle in the form of barrier crossing has been considered by e.g. Beck et al.³⁰ Doye et al.³¹ observed barriers between solid and molten phases for LJ₅₅, a cluster of 55 Lennard-Jones interacting particles. This problem is even more evident in bulk simulations with periodic boundary constraints. Belonoshko and coworkers found that in simulations of bulk aluminum oxide, the melting point was overestimated by several hundred kelvin due to the barrier between the phases.³² Two-phase simulations^32,33 were developed in order to reduce the barrier problem for bulk systems. In such simulations, a cube of solid structure is in contact with a cube of a liquid structure. The interface is a seed for both melting and freezing and hopefully the system will transform into the correct phase for all temperatures. Such a method is however not useful in simulations of small clusters, which should switch phase constantly and not only once. On the other hand, an advantage of small cluster simulations is that the simulations are reliable if frequent phase switches are observed. A corresponding criterion for successful determination of the melting point in bulk simulations is more difficult to define.

What is the character of the phase barrier? We have observed some crossings where the barrier is not a potential energy barrier. When zooming into the potential energy curves in Fig. 10, it is clear that the clusters do not need to overcome a high potential energy barrier. In the solid to molten phase direction, a high step in potential energy must be climbed but in the opposite direction, no potential barrier hinders the transformation. Since the free energy is constituted of a potential energy part and an entropic part, the barrier in the molten to solid direction must be of entropic character only. Zwanzig and Zhou described in a simple model how an exclusively entropic barrier might very well represent the decay behavior of the transformation between different states.⁴⁹ The dynamics of the melting might be illustrated by the time evolution of the distance from the centre of mass in a simulation. In Fig. 12, the distance to the centre of mass for the 55 atoms in Pd₅₅ is drawn versus number of macro steps. The temperature is 900 K which is close to the melting point. To the left, the cluster is in the solid phase and one particle is in the middle, 12 particles are at about 2.6 Å from the centre of mass, and the rest of the atoms are in the outermost shell. At about 400 000 macro steps, one atom in the outermost shell moves out of the shell and becomes a so-called “floater”. (The term was originally used by Cheng et al.⁵⁰) This atom is shown with a thick grey line to distinguish it from the other atoms. This floater causes the first melting but the cluster immediately switches back to the solid phase again. The second melting, at about 600 000 macro steps is not preceded by floaters but is very sudden in character. Thus two kinds of melting is observed. In the simulations, the melting with a preceding floater seems to be most common. All the freezing events (at 450 000, 750 000, 820 000 and 850 000 macro steps) are very sudden. The melting seems to be triggered by floaters or defects but we have not observed any special event just before freezing.

Fig. 12 The distance to the centre of mass for each of the 55 atoms in Pd₅₅ in the course of a simulation is shown at 900 K. One of the atoms is indicated by a thick grey line. (The figure corresponds to Fig. 10).

The barrier between the phases will cause hysteresis of the caloric curve. In the simulations presented this far, the initial configurations were always solid and the barrier hinders the solid to molten phase switch resulting in too high an estimated melting point. In contrast to the MC-CAN simulations, we have seen molten to near-solid transformation for Pd₁₄₇ in MC-ANNEAL simulations. The cluster was annealed from T_initial = 1400 K to T_final = 300 K during 100 million steps. The two caloric curves exhibiting hysteresis are shown in Fig. 13. As a matter of fact, the molten cluster never reaches the true solid configuration but a glass configuration of a little higher potential energy. The barrier problem seems to increase with cluster size. The barrier between the phases is indeed expected to increase with cluster size.⁵¹

Fig. 13 Total energy per atom is plotted versusT for Pd₁₄₇. The solid line with the triangles is the result from the MC-CAN simulation. The grey line is from an MC-ANNEAL simulation where the temperature was reduced by 11 × 10⁻⁶ K per step.

A crucial question is how to solve the problem with barriers in simulations? In the MC simulation, the maximal distance an atom can move in one step is set so that 50% of all attempts are rejected. Is it possible that it is easier for the cluster to overcome the phase barrier if the atoms move by shorter or longer steps? In order to test this we simulated Pd₅₅ at 900 K with maximal displacements giving 35% and 65% rejection, respectively. The simulations did not indicate that another maximal displacement would enhance the number of phase switches. In the literature, different techniques are proposed to ease the overcoming of barriers. In the j-walking technique used by e.g. Moseler and Nordiek⁸ and Frantz,¹² configurations corresponding to a higher temperature are occasionally incorporated in the simulation. With a modification of this method called q-jumping MC, Calvo et al.¹³ found that a solid-molten phase barrier was indeed easier to overcome. In a future publication we will address the problem with phase barriers in clusters with another refined simulation method: the reference system equilibration method.⁵²

VI. Results for Pd₃₄ and Pd₅₄

This far we have only considered icosahedral clusters with closed geometric shells. In this work we have also studied what happens to the melting point when one atom and half the outer shell, respectively, are removed from Pd₅₅. We expect the melting points to be lower than for Pd₅₅, since when the outer shell is not closed, the mobility of the atoms is increased. The thermal behavior of the Pd₅₄ cluster is very similar to that of Pd₅₅. The caloric equation was fitted to the model DOS, giving the parameters in Table 1 and Table 2. As expected, the melting point decreases slightly compared to Pd₅₅.

In contrast, the thermal behavior for Pd₃₄ is different. In Fig. 14, the energy, C_v, δ and the radii of gyration from MC-CAN simulations are drawn versus the temperature and in Fig. 15 the caloric curves from MC-MICRO simulations are drawn. It is remarkable that there is no peak in the heat capacity and for this small cluster the reason cannot be a high barrier between phases. The heat capacity does however show two regions: one for a solid phase and one for a molten phase. The lack of peak in-between, is interpreted as indication that there is no coexistence between phases with significantly different energies. The difference between the solid and molten phases for Pd₃₄ is very small. That is supported by the lowest panel in Fig. 14, where the energy does not show the typical S-shape. The increasing slope only indicates that the atoms probe more and more anharmonic parts of the potential energy surface. An evolution of the kinetic energy without the S-bend is also found in Fig. 15. Furthermore, in the upper panel of Fig. 14, the radii of gyration do not indicate any qualitative change of shape upon melting. From the plot of δ we can see that the atoms start to switch position at about 350 K (MC-CAN) and E = −2.76 eV (MC-MICRO, not shown in figure). For Pd₃₄ we regard the onset of atomic mobility as the melting. An equivalent melting phenomenon was reported for Al₄₃ by Sun et al.⁵³

Fig. 14 Results are shown for Pd₃₄versusT. From above: the maximal, medium and minimal radii of gyration, δ, C_v and the total energy per atom.

Fig. 15 The caloric curves of the Mc-MICRO simulations are shown for Pd₃₄. The rings show the kinetic energy and the pentagrams show the potential energy minus E₀versus total energy.

VII. Conclusions

We have performed Monte Carlo simulations of equilibrating palladium clusters in the canonical ensemble and microcanonical ensembles. In order to avoid the natural fragmentation of a single cluster in an infinite volume at fixed temperature, we have used a new definition of a non-fragmented cluster; the cluster has not fragmented if all atoms are connected via chains of adjacent atoms. Although fragmentation in principle must be prevented for all temperatures, no palladium cluster actually fragmented in reasonable time below 1500 K. The clusters stayed intact sufficiently long to study their equilibrium properties.

By use of a two-state model for the density of states, the simulation results from the microcanonical simulations were converted into temperature-dependent data also for the smallest clusters. The agreement between canonical and microcanonical simulations as well as molecular dynamics simulations with Nosé-Hoover thermostat by Grönbeck et al.²² proved to be excellent with a few exceptions in the coexistence region. In the melting region, the accuracy of δ seems to be worse than for other quantities. Therefore we recommend quantities as the potential energy to be used instead of δ for the calculation of the melting point.

Disregarding that properties of clusters may change dramatically when adding one atom, a rule of thumb quoted by many authors is that the melting point increases with size and finally approaches the bulk melting point.²⁸ In our simulations, Pd₁₃ is an exception from this rule, having the highest melting point of the clusters in this study. However, in support of our findings we note that Kusche in recent experiments found that Na₅₇ has a higher melting point than all the cluster sizes they studied in the size interval from 60 to 200 atoms.¹⁸

The problems associated with free energy barriers between the solid and molten phases that has previously been reported^30,31 were observed in the simulations of Pd₅₄, Pd₅₅, Pd₁₄₇ and Pd₃₀₉. The phase barrier hinders correct calculation of the equilibrium constant in the coexistence region and thus the melting point has not been accurately determined for these clusters. With annealing simulations we can however find a lower limit to the melting point, but as in the case of Pd₁₄₇, with a reasonable cooling rate the width of the hysteresis is too large to give a useful estimation.

If a simulation can be run so that frequent phase switches occur, the phase barrier does not hinder ergodicity and the behavior within the coexistence region is reliable. In the case of Pd₃₀₉, no such fluctuation was ever observed and the melting of this cluster may not be even qualitatively described. The cluster may behave like the smaller clusters without mixing of the two phases as the simulations at a first glance suggest. However, the cluster may also be so large that in the microcanonical ensemble it should, in principle, melt gradually with increasing total energy but in order to overcome the barrier, the total energy must be so high that the equilibrium state is the totally molten cluster. Half-molten configurations may thus be artificially hidden by the barrier. We have studied the crossing of the phase barrier of Pd₅₅ in more detail. At melting a step upwards in potential energy must be taken, while at freezing the barrier is of entropic character only.

In order to get results that could successfully be compared with experimental values the interaction between the atoms must be more accurately modeled. However, of the same significance is to solve the problem with phase barriers. Methods to more efficiently overcome barriers in simulations have been developed and used, for instance the j-walking⁸ and the q-jumping variations of the MC method.¹³ These methods are useful and if they satisfy the criterion that phase fluctuations are observed, the results are reliable. Methods aiming to determine the melting point by computing the DOS of the isolated isomers have been developed.^4,3,10 Although the estimated number of molten subisomers is enormous (e.g., 10²¹ for LJ₅₅⁵⁴), these methods are promising. In a future publication, we will return to one of the methods available to directly calculate the DOS.

Finally we note that the melting of Pd₃₄ is of different character than that of the other clusters. No distinction is observed between the solid and molten structures and we consider melting to occur when the atoms become mobile in the cluster. The onset of mobility is in this case at a considerably lower temperature than for the other clusters.

Appendix: Definition of a fragmented cluster

We have used the following definition of a non-fragmented cluster: A configuration of atoms form a cluster if all atoms are linked to all other atoms via a chain of adjacent atoms. Two atoms which are closer to each other than an element-specific limit B are said to be adjacent.

Let us define the matrix A so that

where r_ij is the distance between atoms i and j. Then calculate the matrix M = A^N−1, where the matrix multiplication is either standard or Boolean. The following definition of a cluster is then equivalent to the one above: If all entries in M are non-zero, the configuration is a cluster, but if any element is zero it is not. The proof is omitted but the result follows physically from the fact that a matrix element in M can only be zero if there is no way to walk between the corresponding atoms in less than N steps between neighboring atoms. The matrix calculation is simple but somewhat time-consuming. In the simulations however, we have the extra information that before the move of one atom, the atoms did form a cluster. With this knowledge, the determination of fragmentation can be done more efficiently.

In order to estimate a reasonable value of B, the probability distribution of configurations for Pd₂ and Pd₃ were used. These probability distributions show minima when the interatomic distances are 4–5 Å. The configurations of Pd₂ and Pd₃ with interatomic distances less than 4 Å would therefore be of non-fragmented character and when any distance is greater than 5 Å, atoms have evaporated. With this indication of a reasonable B, simulations of Pd₁₃ were performed for various B between 4.0 Å and 5.8 Å. The simulations for different B gave almost identical results and all the simulations can safely be considered to represent “non-fragmented” clusters. The results in Secs. 5 and 6 are all for B = 4.0 Å.

Acknowledgements

We are grateful for assistance from Anatoly Belonoshko (Uppsala University) and Henrik Grönbeck (Göteborg University). The project was financially supported by The Swedish Foundation for Strategic Research, SSF and The Swedish Research Council, NFR.

References

1. R. S. Berry, T. L. Beck, H. L. Davis and J. Jellinek, Adv. Chem. Phys., 1988, 52, 75.
2. H.-P. Cheng, X. Li, R. L. Whetten and R. S. Berry, Phys. Rev. A, 1992, 46, 791.
3. P. Labastie and R. L. Whetten, Phys. Rev. Lett., 1990, 65, 1567.
4. D. J. Wales, Mol. Phys., 1993, 78, 151.
5. T. L. Beck, J. Jellinek and D. J. Wales, J. Chem. Phys., 1987, 87, 545.
6. M. Bixon and J. Jortner, J. Chem. Phys., 1989, 91, 1631.
7. R. E. Kunz and R. S. Berry, Phys. Rev. Lett., 1993, 71, 3987.
8. M. Moseler and J. Nordiek, Phys. Rev. B, 1999, 60, 11 734.
9. J. P. Neirotti, F. Calvo, D. L. Freeman and J. D. Doll, J. Chem. Phys., 2000, 112, 10 340.
10. L. M. Amon and W. P. Reinhardt, J. Chem. Phys., 2000, 113, 3573.
11. A. Rytkönen, S. Valkealahti and M. Manninen, J. Chem. Phys., 1998, 108, 5826.
12. D. D. Frantz, J. Chem. Phys., 1995, 102, 3747.
13. F. Calvo and F. Spiegelmann, J. Chem. Phys., 1999, 112, 2888.
14. J. Jellinek, T. L. Beck and R. S. Berry, J. Chem. Phys., 1986, 84, 2783.
15. C. L. Cleveland, W. D. Luedtke and U. Landmann, Phys. Rev. Lett., 1998, 81, 2036.
16. M. Schmidt, R. Kusche, W. Kronmüller, B. von Issendorff and H. Haberland, Phys. Rev. Lett., 1997, 79, 99.
17. M. Schmidt, R. Kusche, B. von Issendorff and H. Haberland, Nature, 1998, 393, 238.
18. R. Kusche, Th. Hippler, M. Schmidt, B. von Issendorff and H. Haberland, Eur. Phys. J. D, 1999, 9, 1.
19. J. Westergren, H. Grönbeck, S. G. Kim and D. Tománek, J. Chem. Phys., 1997, 107, 3071.
20. J. Westergren, H. Grönbeck, A. Rosén and S. Nordholm, J. Chem. Phys., 1998, 109, 9848.
21. J. Westergren, H. Grönbeck, A. Rosén and S. Nordholm, Nanostruct. Mater., 1999, 12, 281.
22. H. Grönbeck, D. Tománek, S. G. Kim and A. Rosén, Chem. Phys. Lett., 1997, 264, 39.
23. X. Chen, J. Zhao, Q. Sun, F. Liu, G. Wang and X. Shen, Phys. Status Solidi B, 1996, 193, 355.
24. L. J. Gallego, M. J. Grimson and C. Rey, Phys. Rev. B, 1995, 51, 5518.
25. M. J. Grimson, Chem. Phys. Lett., 1992, 195, 92.
26. J. Akola and M. Manninen, Phys. Rev. B, 2001, 63, 193 410.
27. F. Ercolessi, W. Andreoni and E. Tosatti, Phys. Rev. Lett., 1991, 66, 911.
28. A. Rytkönen, S. Valkealahti and M. Manninen, J. Chem. Phys., 1997, 106, 1888.
29. S. Nosé, J. Chem. Phys., 1984, 81, 511; W. G. Hoover, Phys. Rev. A, 1985, 31, 1695.
30. T. L. Beck and R. S. Berry, J. Chem. Phys., 1988, 88, 3910.
31. J. P. K. Doye and D. J. Wales, J. Chem. Phys., 1995, 102, 9673.
32. A. B. Belonoshko, Geochim. Cosmochim. Acta, 1994, 58, 4039; A. B. Belonoshko and L. S. Dubrovinsky, Am. Mineral, 1996, 81, 303; R. Ahuja, A. B. Belonoshko and B. Johansson, Phys. Rev. E, 1998, 57, 1673; A. B. Belonoshko, Phys. Chem. Minerals, 1998, 25, 138.
33. J. D. Kubicki and A. C. Lasaga, Am. J. Sci., 1992, 292, 153.
34. M. Andersson, J. L. Persson and A. Rosén, J. Phys. Chem., 1996, 100, 12 222.
35. H. Haberland, Clusters of Atoms and Molecules I, Chemical Physics, Springer-Verlag, Berlin, 1995, 52.
36. G. Scoles, Atomic and Molecular Beam Methods, Oxford University Press, Oxford, 1992, vol. 2, p. 231.
37. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, J. Chem. Phys., 1953, 21, 1087.
38. J. K. Lee, J. A. Barker and F. F. Abraham, J. Chem. Phys., 1973, 58, 3166.
39. S. M. Thompson, K. E. Gubbins, J. P. R. B. Walton, R. A. R. Chantry and J. S. Rowlinson, J. Chem. Phys., 1984, 81, 530.
40. L. Råde and B. Westergren, Beta-Mathematics Handbook, Studentlitteratur, Lund, 3rd edn., 1995.
41. E. S. Severin, B. C. Freasier, N. D. Hamer, D. L. Jolly and S. Nordholm, Chem. Phys. Lett., 1978, 57, 117.
42. H. W. Schranz, S. Nordholm and G. Nyman, J. Chem. Phys., 1990, 94, 1487.
43. D. Tománek, S. Mukherjee and K. H. Bennemann, Phys. Rev. B, 1983, 28, 665; W. Zhong, Y. S. Li and D. Tománek, Phys. Rev. B, 1991, 44, 13 053.
44. G. L. Estiù and M. C. Zerner, J. Phys. Chem., 1994, 98, 4793.
45. C. Rey, L. J. Gallego, J. García-Rodeja, J. A. Alonso and M. P. Iñiguez, Phys. Rev. B, 1993, 48, 8253.
46. D. Frenkel and B. Smit, Understanding Molecular Simulation, Academic Press, New York, 1996.
47. F. Calvo and P. Labastie, Chem. Phys. Lett., 1995, 247, 395.
48. D. J. Shaw, Introduction to Colloid and Surface Chemistry, Butterworth-Heinemann, Oxford, 1992, p. 64.
49. H.-X. Zhou and R. Zwanzig, J. Chem. Phys., 1991, 94, 6147; R. Zwanzig, J. Chem. Phys., 1992, 97, 3587.
50. H.-P. Cheng and R. S. Berry, Phys. Rev. A, 1992, 45, 7969.
51. R. M. Lynden-Bell and D. J. Wales, J. Chem. Phys., 1994, 101, 1460.
52. L. Ming, S. Nordholm and H. Schranz, Chem. Phys. Lett., 1996, 248, 228.
53. D. Y. Sun and X. G. Gong, Phys. Rev. B, 1998, 57, 4730.
54. J. P. K. Doye and D. J. Wales, J. Chem. Phys., 1995, 102, 9659.

---

Footnote

† Present address: AstraZeneca R&D Mölndal, SE-431 83, Mölndal, Sweden.

---