Melting dynamics and isomer distributions of small metal clusters

The structure and dynamics of NiN, AgN and AuN (N = 6–30) clusters have been studied extensively by Monte Carlo and molecular dynamics methods based on the Sutton–Chen many-body potential. An exhaustive search for low-energy minima on the potential energy surface was carried out using the eigenvector-following technique. The exponential increase in the number of isomers with atomic size is demonstrated and compared. The binding energies and point groups of global minimum and first two isomers of NiN, AgN and AuN (N = 6–14) clusters are listed. The melting properties and temperatures of the clusters are reported.


Introduction
Clusters are of fundamental interest both due to their own intrinsic properties and because of the central position that they occupy in molecular and condensed matter science. Clusters are distinct from bulk materials because they are dominated by surface species and, consequently, clusters have structures and properties that often differ from anything that can be observed in bulk [1]- [3]. Structural and dynamical properties of atomic and molecular clusters are the subject of increasing intense research activity both theoretically [4]- [24] and experimentally [25]- [36]. Small metal clusters are important in homogeneous nucleation [37,38], crystal growth [1,39], catalysts [40]- [42], understanding the structure of amorphous materials. They are also important due to their potential applications in nanoelectronic and memory devices [43]- [45].
When the cluster structure is considered, it is important to realize that clusters exist within the domain of a complex many-particle potential surface with many local minima. Each of the minima can be associated with a cluster structure that can be thought of as an isomer of the potential energy function. Many methods have been developed to explore the potential surfaces of clusters including the local minima of the surface and the transition states. Some of these methods have been simulated by annealing approaches [46,47], conjugate-gradient techniques [48,49], genetic algorithm [50]- [54] and eigenvector-following methods [55]. There are also ab initio calculations to find out the global minimum and low-energy isomers of small metal clusters [56]- [61]. These first-principles calculations for understanding the structure and electronic properties of clusters are especially focused on small gold clusters and their anions. Landman and co-workers performed DFT calculations, supported by the photoelectron spectroscopy experiments and found two-dimensional structures for N < 13 gold anions and proposed two-dimensional structures for the neutral clusters of the same size [56,60,61].
In recent years, there has been an increased interest in investigating the physical properties of metallic clusters as a function of energy or temperature. These canonical and microcanonical simulations have shown that some properties of clusters as a function of temperature or energy exhibit rapid changes that correspond to changes seen in bulk systems undergoing phase transitions. It has become common to identify these regions of temperature as cluster-phase transition analogues [62]. Most of the work on phase changes in clusters has focused on analogues of melting transitions. A comprehensive understanding of the melting phenomenon is crucial because it is directly related to the stability, structure and size of the cluster. Melting in a macroscopic object occurs at some well-defined temperatures. This is no longer true for a small particle or cluster. It is well known that the melting point of a small cluster with free surface, i.e., an isolated cluster in vacuum, is lower than that of the infinite bulk solid [63]. This can be understood according to phenomenological models by emphasizing the enhanced surface effects due to the small size.
Recently, Haberland and co-workers obtained caloric curves from the photo-dissociation measurements of sodium clusters and showed that the melting point and latent heat of fusion vary non-monotonically with size [31,32,36]. Even though some theoretical works have been carried out to explain the experimental results, a full explanation has not been achieved up to now [64]- [69].
In this paper, we have studied isomer distributions of Ni N , Ag N and Au N (N = 6-14) clusters, and the ground state structures and the melting dynamics of Ni N , Ag N and Au N (N = 6-30) cluster by Monte Carlo and molecular dynamics methods. In the following section, the Sutton-Chen potential and the details of the numerical procedures are presented. In section 3, we discuss the results. Finally, we present a brief summary of our study.

The potential and computational procedure
Following Finnis and Sinclair [70], the potentials have the form where the first term represents the repulsion between atomic cores and the second term models the bonding energy due to the electrons. As the electrons are not included explicitly, the local density, ρ i , in the second term is the local density of atoms. In Sutton and Chen's [71] potentials, the local density in the second term is given by and the pair repulsion is also modelled by a reciprocal power so that the complete potential is Although this is truly a many-body potential, the force on each atom can be written as a sum of pairwise contributions where Here r ij is the separation between atoms i and j, C is a positive dimensionless parameter, ε is a parameter with the dimensions of length and m and n are positive integers. The three parameters ε, C and a are not independent and are completely determined by the equilibrium lattice parameters and lattice energy of the face-centred lattice.
In our calculations, we use the same parameters as given by Sutton and Chen [71] for nickel, silver and gold metals (table 1). This potential has been shown to produce bulk and surface properties quite accurately as reported in [5,72,73] and references therein.
The groundstate geometrical structures of Ni N , Ag N and Au N (N = 6-30) clusters and isomer distributions of Ni N , Ag N and Au N (N = 6-14) clusters are studied by eigenvectorfollowing [55,74] technique using Monte Carlo simulations. The starting configurations of each cluster was achieved by placing the atoms randomly in a box of size 10 × 10 × 10 Å. The distances between atoms were prevented from being too small or too large. To obtain all the isomers of clusters that can be predicted by this technique, 3 × 10 3 distinct random configuration sets for each cluster size of each element were prepared. Then these sets of random configurations were optimized by eigenvector-following technique [55]. This procedure is repeated for the new sets of distinct random configurations and it is stopped when the total number of isomers of each cluster remains constant. During this process, it is observed that for the atom numbers N = 6-9, 9 × For N = 6-9, all three metals have the same three (two for N = 6) lowest-energy isomers and in the same order of energy. For N = 6, we found only two lowest-energy isomers with the Sutton-Chen potential but three lowest-energy isomers were found with the Gupta potential [9]. The global minima and isomers are in agreement with the results obtained by the Gupta potential for N = 6 and 7. The Ni and Ag cluster global minima are based on variant of the icosahedral structure for N = 12-14. For Ni and Ag, at N = 12-14 and for Au at N = 13, 14 the global minima found here with the Sutton-Chen potential are identical to those found with the Gupta potential. The global minimum of Au 12 cluster found with the Gupta potential corresponds to the structure that we found here for the first isomer of Au 12 cluster with the Sutton-Chen potential.
Since there is a lack of information about the point group in [9], it is difficult to compare the isomers of the clusters for N = 12-14. All the global minimum structures of Sutton-Chen Ni N , Ag N andAu N (N = 6-30) clusters which were obtained here by eigenvector-following technique, match with the previously published results for those Sutton-Chen clusters obtained by the basin hopping algorithm [73].
In the solid-liquid phase transition calculations described in this paper, the Hamiltonian equations of motions are solved for the Sutton-Chen potential for all the atoms in a cluster on a grid of total energies using Hamming's modified fourth-order predictor-corrector propagator with a step size of 0.5 × 10 −15 s. The clusters are prepared with zero initial total and angular momentum. Then, the clusters are gradually heated up to 2000 K. Trajectories of length of 5 × 10 5 steps are generated on a grid of total energies that are large enough to observe the solidliquid transitions in the clusters. The total energy in the individual runs is conserved within 0.001%.
To analyse the melting-like transitions of clusters, two different types of magnitudes are used: (1) global quantities such as temperature and the rms bond length fluctuation (δ) (these quantities are calculated as time averages over a whole trajectory at a given energy) and (2) time-dependent quantities such as the mean squared displacement and the velocity autocorrelation function. These quantities can be calculated as averages over trajectories (all of them corresponding to the same energy) for each time, or equivalently, as averages over well-separated time origins defined along a single trajectory.

Global quantities
The cluster temperature is given by, where n is the number of atoms in the cluster, k B is the Boltzmann constant, E k is the total kinetic energy of the cluster and represents the time averages along the whole trajectory.   The relative rms bond-length fluctuation (δ) is defined by This quantity experiences an abrupt (but continuous) increase associated with the melting-like transition of the cluster. As an estimation of the melting temperature T m , the Lindemann criteria (δ > 0.10) is considered.

Time-dependent quantities
The mobility of the atoms in the cluster can be analysed in terms of the mean-squared displacements, where n t is the number of time origins (t oj ) considered along a trajectory. When the atoms have an oscillatory motion around their equilibrium positions, this quantity is a constant as a function of time. As the energy of the cluster increases, the atoms begin to move in a diffusive way and then, the mean-squared displacement increases linearly with time.
The velocity autocorrelation function given by provides complementary information on the motion of the atoms in the cluster. For low energies, the motion of the atoms is highly correlated and C(t) exhibits oscillations as a function of time. This is the typical behaviour of a solid-like cluster. For high energies (corresponding to liquid-like clusters), the oscillations in C(t) are completely lost (after the first minimum) which indicates an uncorrelated motion of the atoms in the cluster. The power spectrum, I(ω), which is the Fourier transform of C(t) and is expressed as where ω is the cyclic frequency and C(t) is the velocity autocorrelation function of the system.

Results and discussion
The ground state structures and isomer distributions of the Ni N , Ag N and Au N (N = 6-14) clusters are obtained and the variation in the number of isomers as a function of atomic size is plotted in figure 4. It is observed that, as the number of atoms increases linearly, the corresponding  number of isomers increase exponentially. The exponential increase in the number of isomers was pointed out before by Stillinger and Weber [75,76]. The graphs in figure 4 are fitted to exponential growth and the equation,   figure 4, gold clusters have more isomers than nickel and silver clusters and also the difference in energy between the global minimum and low-lying isomers are very small for gold clusters with respect to nickel and silver clusters as shown for 13 atomic clusters in figure 5. This situation is the same for other clusters presented here. The energy density of stable states near the global minimum is a function of the attractive range (m values in table 1) of the potential for these metals. In table 1, the m values of nickel and silver clusters are smaller than the m values of a gold cluster. As it is known, as the value of m increases the potential energy function becomes short range. The shorter the range, the higher the density of states and the smaller the energy gap between the global minimum and the higher energy isomers. If we look at figure 5 carefully, it can be seen that at certain total energy values there is a small jump in energy value for Ni 13 and Ag 13 clusters. However, this situation is not observed for the Au 13 cluster. Because the potential energy function for the gold cluster is of short range and so the potential energy surface is more complex, the gold clusters have more density of stable states and have more isomers than the nickel and silver clusters. Therefore, the energy differences between the gold isomers are small and there is no jump in energy value for the Au 13 cluster in figure 5.
By monitoring the melting process, the short time averages of kinetic energies are plotted and from these curves, the potential wells visited by the cluster are examined by quenching all the configurations corresponding to the energy values just before the phase transitions. In this way, isomers are obtained. It is observed that, during the melting process Ni N , Ag N and Au N (N = 6-9) clusters spend more time in their first isomer than the others do. However, the clusters with N > 9 atoms spend more time in the second, third, etc isomers rather than in the first isomers before melting. To illustrate the situation for 13 atomic clusters, it is found commonly that clusters spend much more time in their fourth isomer for Ni 13 , fifth isomer for Ag 13 (as shown in figure 8) and seventh isomer for Au 13 . Now, to discuss the melting behaviour of the clusters, firstly we look at the variation of the total energy of clusters as a function of temperature. There are two main types of caloric curve observed for clusters in the size range of N = 6-30 atoms, as shown in figure 6 for Ag 13 and Ag 14 respectively. For the clusters with closed shells such as Ag 13 , the change in slope at a characteristic temperature is much more abrupt than in the open-shell clusters such as Ag 14 . This can be seen more clearly in figure 7, where we have plotted the rms bond-length fluctuation (δ) with respect to temperatures.
For the closed-shell cluster (i.e. clusters with geometries where the outer shells are fully decorated with atoms) Ag 13 , the δ versus T plots show three distinct regions. In the first region that corresponds to solid phase, δ varies slowly with temperature until it shows a sharp change. In the second region, δ varies rapidly with temperature and melting starts in the system. The third region corresponds to the liquid state where again δ varies slowly with temperature.
For the open-shell clusters, the variation of δ versus T suggests that these clusters melt in two stages. In the first stage, δ shows an abrupt change at a rather low temperature (of about 400 K for Ag 14 ). In the second stage, at a higher temperature (∼ 800 K for Ag 14 ), δ shows again a jump in its value over a relatively broad range of temperatures. This feature of two-stage melting is seen to be a characteristic of all open-shell nickel, silver and gold clusters such as N = 14-18 and 20-30. The pre-melting behaviour is strongly influenced by vacancies and anharmonicities that are distinctly larger in open-shell clusters than in closed-shell clusters. The surface atoms, which do not belong to the complete shell, are in general, more weakly bound than the atoms in the cluster, and, therefore, the amplitude of these atoms are significantly larger than those in closed-shell clusters. This gives rise to enhanced anharmonicities, and provides fluidity in the In the lower panel, the cluster structure corresponds to global minimum (I h ) and frequently visited isomers (C 1 ) just before phase transition.
system at lower temperatures, a fact that leads to pre-melting behaviour in these systems. The pre-melting behaviour is very similar to the surface melting phenomena observed in bulk crystals [77] and arises due to the difference in coordination numbers between the atoms on the surface and the atoms in the bulk. The solid and liquid characters are further confirmed by calculating the short time averages of kinetic energy, velocity autocorrelations and diffusion coefficients. Figures 8 and 9 display the instantaneous temperatures, i.e. short time-averaged temperatures as a function of time for different total energies per atom for Ag 13 and Ag 14 . The time dependence of the instantaneous temperature of a cluster is an informative characteristic of the state of this cluster at a given total energy. As the total energy increases, the pattern of the graphs changes from constant to moderately fluctuating and then to highly fluctuating. The graphs (a) correspond to a lowenergy rigid cluster, the graphs (b) to clusters that are undergoing structural changes and (c) to a high-energy liquid-like cluster. In these figures, graph (b) of Ag 13 shows wide potential wells that suggest that, at this stage, the cluster undergoes structural changes between the energetically close isomers. For the Ag 14 cluster, the graphs (a) show a unimodal distribution, whereas graph (b) has a bimodal character. Graph (b) has narrow but sharp peaks, i.e., the cluster visits energetically close isomers for very short time intervals. This stage resembles the fluctuating state described by Sawada and Sugano [78,79] because in this state there is no pronounced diffusion as shown in figure 10; this type of behaviour of clusters is called pre-melting. Graph (d) corresponds to high-energy liquid-like clusters.
In figure 10, to corroborate the melting-type nature of the transition, the mean-squared displacements of the atoms in the cluster are displayed respectively for two clusters (Ag 13 for Ag 14 , the velocity autocorrelation functions have significant contribution from very low frequency modes due to their liquid character that is shown clearly in the Fourier transforms. The overall change in the pattern of Fourier transform graph with energy (broadening and shift towards lower frequencies) is typical of a solid-to-liquid-like transition [80,81].  The melting temperatures of the Ni N ,Ag N andAu N (N = 6-30) clusters have been estimated by monitoring the caloric curve, rms bond-length fluctuation, diffusion coefficients and velocity autocorrelation functions and their Fourier transforms. Generally, the Lindemann criterion is considered in determining the melting temperatures of the clusters. Similar criteria for melting in small clusters have also been used previously [22,23]. The size dependence of the melting temperatures of the Ni N , Ag N and Au N (N = 6-30) clusters are shown in figure 13. For the clusters showing pre-melting, the pre-melting temperature is chosen as a measure of the melting point. In these cases, while it is relatively easy to determine the temperature at which the melting starts in the clusters in the first stage (first jump in δ), the second stage (second jump in δ) of melting proceeds over a relatively wider range of temperature, thereby making the determination of a second set of melting temperatures quite difficult. When the melting temperatures of the three metals are compared, it can be concluded that the Ag and Ni clusters have similar melting behaviour. It is also observed that the clusters with N < 13 atoms, in general, have higher melting points than those consisting of more than 13 atoms (except Ni 19 and Ag 19 ). For the Ni 13 , Ag 13 , Ni 19 and Ag 19 clusters, there is a sharp increase in melting temperature, since those clusters have highly stable structures compared to other sizes. However, even the Au 13 cluster has an icosahedral structure and there is no sharp increase in melting temperature as for the Ni 13 and Ag 13 clusters. This may be explained by the existence of the small energy difference between the global structure and the next higher-energy structure as indicated in table 2. Similarly, the Au 19 cluster has glass-like structure with only one internal atom and is not stable as the Ni 19 and Ag 19 clusters. It is interesting to note that for clusters with N < 13, all the atoms are surface atoms, while in clusters N 13, there is at least one atom that can be identified as a bulk atom and they have open-shell structures. Therefore, the atom in the open shell moves easily and the cluster shows premelting behaviour at the lower temperatures.

Summary
The global minimum structures and isomer distribution of clusters have been studied by the eigenvector-following technique using Monte Carlo simulations, and the binding energies and point groups of global minima and first two isomers of the Ni N , Ag N and Au N (N = 6-14) clusters are listed. It is found that the gold clusters have more isomers than nickel and silver clusters and the energy differences between the global minimum and the other isomers are very small for gold clusters with respect to nickel and silver clusters. This is because of the higher density of states of gold clusters, and the attractive part of the potential energy function for gold is short-ranged than nickel and silver. The number of isomers for clusters presented here increases exponentially with cluster size.
Molecular dynamics simulations are employed using Sutton-Chen potential for the Ni N , Ag N and Au N (N = 6-30) clusters and melting temperatures of clusters are obtained. It is shown that small clusters undergo a transition from a solid to liquid-like state with an increase in internal energy. The frequently visited isomers of the Ni N , Ag N and Au N (N = 6-30) clusters are obtained by quenching all the configurations corresponding to the energy values just before the phase transitions. The melting temperatures are found to be a non-monotonic function of the cluster size, with some clusters exhibiting 'pre-melting' behaviour in analogy with surface melting in bulk crystals. In conclusion, the melting temperatures of clusters which have N < 13 atoms and closed-shell clusters N = 13, 19 (for nickel and silver) are higher than those for other clusters.