- Email: [email protected]
Molecular dynamics simulation of the melting-like transition in K1Na54
Computational Materials Science 35 (2006) 174–178 www.elsevier.com/locate/commatsci
Molecular dynamics simulation of the melting-like transition in K1Na54 Andre´s Aguado
*,1,
Jose´ M. Lo´pez, Sara Nu´n˜ez
Department of Theoretical Physics, University of Valladolid, Valladolid 47011, Spain Received 24 April 2004; accepted 11 October 2004
Abstract A set of molecular dynamics simulations at constant total energy has been performed in order to investigate the melting-like transition of the impurity-doped cluster K1Na54. An orbital-free density-functional-theory technique is employed to obtain the forces on atoms in an efficient way. The total simulation time is approximately 1.5 ns, which is required to obtain statistically meaningful thermal averages. The presence of just one K impurity reduces considerably the melting temperature value (as compared to Na55). Although only one peak is observed in the thermal evolution of the specific heat, analysis of the diffusion constants shows that the cluster melts in two separate steps: Na atoms are able to diffuse at lower temperatures than the K impurity. 2005 Elsevier B.V. All rights reserved. PACS: 31.15.Qg; 36.40.Ei; 36.40.Sx; 64.70.Dv Keywords: Phase transitions in clusters; Alkali clusters; Ab initio molecular dynamics
1. Introduction A series of recent experiments by Haberland and coworkers [1a–d] has renewed the interest in the melting transition of small atomic clusters. These experiments reveal that variation of melting temperature Tm with cluster size N is non-monotonic, contrary to expectations based on classical arguments stating that melting temperature of a finite liquid drop should monotonically decrease as its radius shrinks. The detailed evolution of Tm with N shows enhanced stability of the solid phase for selected ranges of cluster sizes which are close to either electronic and/or geometrical shell closings. Thus, it is now accepted that both electronic and structural
*
Corresponding author. Tel.: +34 83 423147; fax: +34 83 423013. E-mail address: [email protected] (A. Aguado). 1 Supported by the ‘‘Ramo´n y Cajal’’ program (Spanish Ministry of Science and Technology). 0927-0256/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2004.10.011
effects will play a role in explaining the observed results, but apart from this tacit assumption, no definite explanation of the experimental observations has emerged yet. Computer simulation of the melting process may provide physical insight into the microscopic details of the melting-like transition and, hopefully, reproduce and explain the experimental Tm(N) curve. A very recent attempt in this direction has been given by Manninen et al. [2] These authors perform molecular dynamics (MD) simulations of melting in NaN (N = 40–355) and reproduce non-monotonic Tm(N) behavior. However, due to computational limitations, interatomic forces in all clusters with N P 93 are calculated employing a phenomenological potential, and the Tm(N) curve is not in good agreement with the experimental one for N P 142. The reasons for this discrepancy may be: (a) purely electronic effects (not well described by any phenomenological potential) have a strong influence; (b) interatomic forces provided by that potential might be not accurate
A. Aguado et al. / Computational Materials Science 35 (2006) 174–178
enough so that correct ground state (GS) isomers (for example, those predicted by an ab initio methodology) are not reproduced; (c) even if the potential correctly reproduces the ab initio energetic ordering of isomers, it is very difficult to locate the GS isomer when so many atoms are involved; etc. Agreement with experiment is much better for smaller clusters. Specifically, Manninen et al. reproduce the experimental result T m ðNaþ 55 Þ > T m ðNaþ 93 Þ, and MD simulations by Aguado et al. [3a,b] reproduce experimental melting temperatures of selected clusters in the range N = 8–142 to within a few percent. In the future, employment of more accurate interatomic forces and a thorough search for GS isomers may lead to satisfactory results also for large sizes. Consideration of the melting-like transition in impurity-doped metal clusters (and in multicomponent clusters in general) is even more complicate from a computational point of view. First of all, compositional ordering increases significantly the number of different isomer configurations, and thus the search for the GS isomer is much more time consuming. Also, surface segregation of the component with lowest surface energy may result in specific features such as partial melting of one component or incomplete mixing. Finally, and most important for the melting problem, the presence of impurities in an otherwise homogeneous cluster might alter structural properties, influencing thermal properties and, in particular, the location of the melting point. Although most of the work on multicomponent clusters has been devoted to structural properties, a few molecular dynamics simulations exist studying their freezing [4] and melting [5–8] transitions, and the influence of temperature on segregation [9]. In the present paper we report the results of MD simulations of the melting-like transition of the impuritydoped cluster K1Na54. To better identify the details of the solid-to-liquid transition, we have employed an orbital-free (OF) version [10,3a] of density-functional-theory (DFT) [11] (with core electrons replaced by pseudopotentials) in the calculation of interatomic forces, which are then used in a classical (Newtonian) description of atomic motions. Computational expense scales linearly with system size in OF approaches, allowing consideration of MD simulation runs long enough to extract meaningful statistical averages. The necessary employment of local pseudopotentials and approximate density functionals for the electronic kinetic energy reduces the accuracy of these methods, as compared to traditional Kohn–Sham (KS) approaches [12]. Nevertheless, OFMD includes an approximate explicit description of the electronic density, and thus provides a more transferable description of metal clusters than phenomenological potentials, for example. At the same time, a better statistical sampling of phase space (as compared to traditional KS approaches) may be achieved with only moderate effort. We will show that introduction of a
175
single K impurity in a Na55 host significantly affects thermal behavior: (1) the melting temperature is considerably lower than that of homogeneous Na55; (2) melting in K1Na54 proceeds in two separate steps, contrary to Na55, which melts in one single stage [3b]. We will also show the necessity of evaluating several melting indicators in order to obtain a complete picture of the melting process. In Section 2 we briefly present some technical details of the method. The results are presented and discussed in Section 3 and, finally, Section 4 summarizes our main conclusions.
2. Details of computer simulations Our implementation of OFMD is based on an expansion of the electron density in a basis set formed by plane waves. In order to make this implementation efficient, local pseudopotentials are employed to represent the effect of core electrons on the valence density. Thus, only this valence density is explicitly described. Details of the specific pseudopotentials and electronic kinetic and exchange-correlation energy functionals employed are the same as in our previous work [3a,3b]. The cluster is placed in a unit cell of a cubic superlattice, and the set of plane waves periodic in the superlattice is used as a basis set to expand the valence density. ˚ . and an Present calculations use a supercell of edge 35 A energy cut-off in the plane wave expansion of the density of 20 Ry. A 96 · 96 · 96 grid is used to perform Fast Fourier Transforms. The equations of motion are integrated using the velocity Verlet algorithm [13] with a time step equal to 4(3) fs for temperatures lower (higher) than 100 K. These choices resulted in a conservation of the total energy better than 0.1%. The ground state electron density is found for each configuration of atoms visited in the course of an MD trajectory, which allows the accurate calculation of atomic forces by means of the Hellman–Feynman theorem. As a first step, we have tried to locate a plausible GS isomer for K1Na54. Specifically, we have considered complete icosahedral and cuboctahedral structures and all inequivalent positions for the substitutional impurity. In order to explicitly check the reliability of OF predictions, we optimize the same structures employing pseudopotential KS-DFT calculations as implemented in SIESTA [14], applying the same exchange-correlation functional (LDA) for consistency. Standard double-zeta plus polarization (DZP) basis sets were adopted for each atom, together with an energy cutoff of 20 Ry. Stable isomer geometries are found by employing a conjugate gradients optimization technique and requiring the residual force on each atom to be smaller than ˚. 0.001 eV/A Once we have found a candidate GS isomer, we perform a set of constant energy MD simulations in order
A. Aguado et al. / Computational Materials Science 35 (2006) 174–178
Different values of the average temperature were obtained by isokinetic thermalization runs (approximately 5 ps long) previous to each constant energy run. A typical simulation time for the constant energy runs of solid-like clusters is 50 ps, but runs longer than 200 ps were performed for those energies where the melting transition sets in. For energies corresponding to liquidlike clusters, the length of the simulations was always in the range 70–100 ps. The total simulated time for each cluster was approximately 1.5 ns.
3. Results There are four geometrically inequivalent sites in a perfect 55-atom icosahedron. Introducing a substitutional impurity in each of these sites generates four different K1Na54 isomers, whose energies are compared in Table 1. We do not include results for corresponding cuboctahedral isomers as these are located at higher energies. As expected from size mismatch and surface energy arguments, K prefers to be located at a vertex site of the cluster surface. Moreover, OF calculations reproduce the energetic ordering of isomers predicted by SIESTA, which supports the employment of our OF procedure to study thermal properties. From a structural point of view, the only appreciable difference between OF and SIESTA calculations is that OF distances are systematically longer by slightly less than 1%, an error which lies well within the expected accuracy of our methodology. In Fig. 1, we show the caloric curve, specific heat, partial rms bond length fluctuations and diffusion constants of K1Na54, as a function of average temperature. A thermal phase transition is indicated in the caloric
Table 1 Energies (in meV/atom) of several isomers of K1Na54, taking the energy of the corresponding ground state (GS) isomer as a reference
SIESTA OF
a
b
c
d
9.5 10.0
7.6 8.0
GS GS
1.5 1.4
K atom is located at (a) the center of the icosahedron; (b) any of the 12 equivalent positions within the first icosahedral shell; (c) a vertex position of the surface shell (again, the multiplicity of this site is 12); (d) any of the 30 non-vertex sites of the outer shell.
-12.4 -12.45 Na-Na Na-K
-12.5 -12.55 7
0.2 0.15 0.1
δ
0.3 0.25
0.05 0 15
Na K
6
10
5 4
5
D (bohr2/ps)
j
-12.3 -12.35
E (hartree)
to obtain the caloric curve, specific heat, root-meansquared bond length fluctuations and diffusion constants. The way these averages are evaluated is explained in previous works [3a,b]. We just show here explicitly the expression for the ‘‘atomic equivalence indexes’’ [15], as these are less standard quantities X ri ðtÞ ¼ j~ Ri ðtÞ ~ Rj ðtÞj. ð1Þ
Specific Heat (kB)
176
3 0
2 0
100
200
T (K)
300
400
100
200
300
400
T (K)
Fig. 1. Caloric and specific heat curves (left side) and partial rms bond length fluctuations and diffusion constants (right side) of K1Na54, taking the internal cluster temperature as the independent variable.
curve by a change of slope, the slope being the specific heat; the height of the step gives an estimate of the latent heat of fusion. However, melting processes are more easily recognized as peaks in the specific heat or steps in d, as a function of temperature. All three indicators lead to the same value for the estimated melting temperature (130 K), which is a good check for the convergence of our results. Our estimation for the latent heat of fusion is of 0.35 eV. The temperature evolution of the specific heat displays just one main peak centered at Tm = 130 K. At this temperature, there is an abrupt increase in the phase space volume available to the system. No more features are observed in the specific heat, except for temperatures close to 400 K, where evaporation of atoms sets in (incidentally, this last temperature is quite close to the evaporation temperature estimated by Schmidt and Haberland [1d] for Naþ 139 . An examination of the boiling transition is, however, outside the scope of the present communication). Analysis of the temperature evolution of diffusion constants, however, shows that the K atom does not participate in diffusion events until a temperature of 150 K is reached. This is significantly higher than Tm and shows that melting proceeds in two well defined steps: (1) at T = Tm, all sodium atoms (both surface and inner-shell atoms) can freely diffuse across the whole cluster volume, but K atom stays in its vertex position at the surface, without participating in diffusion processes; (2) at T = 150 K, free energy barriers against K diffusion may be overcome and the whole cluster is fluid-like. This two-step melting process is validated by visual inspection of the movement of atoms on a computer screen. The second step is not reflected in the specific heat because the relative increase in the available phase space associated to it is very small. It is not reflected either in the partial bond length fluctuations because diffusion of only Na atoms is enough to induce a step in the thermal evolution of Na–K bond length
A. Aguado et al. / Computational Materials Science 35 (2006) 174–178
fluctuations. It only can be appreciated in indicators which depend on the identity of each single atom, as is the case for the diffusion constants. Fig. 2 shows the time evolution of short-time-averages of the atomic equivalence indexes introduced in previous section, for several representative temperatures. These quantities are adequate to analyze some microscopic details of the melting transition. At T = 120 K, atoms are grouped in shells, representing symmetry equivalent atomic sites. For the case at hand, we find the smallest value of r for the Na atom at the center of the icosahedron; then a group of 12 atoms forming the first icosahedral shell; next, a group of 30 lines representing those atoms of the outer shell capping the edges of the inner Na13 icosahedron; and finally, 12 lines for the atoms at vertex positions of the surface. Within this surface shell, K is clearly identified due to its larger nearest neighbor distances. Although no mixing between different shells is observed, some intrashell disorder already appears at this temperature which is very close to Tm. This disorder is induced by temporary excursions of some atoms to positions of low average radial atomic density. After these temporary excursions, each atom comes back to its original site, that is, this disorder does not involve intrashell diffusion. At a higher temperature (T = 140 K), the amplitude of these motions is large enough to induce mixing between different shells and melting of the whole ‘‘Na sublattice’’. Note that there is not a separate surface-melting stage, as observed for example in Na20 or Na92 [3a,b], but that both surface and inner Na atoms start to diffuse at the same temperature. This is similar to the homogeneous melting process observed for Na55 [3b]. The only difference is that, at this temperature, K atom is not able to diffuse yet. For all temperatures between 140 and 360 K, the atomic equivalence indexes look quite similar, that is, these indicators
σ (a.u.)
1000 900 800
T=120 K
σ (a.u.)
700 600 1000
T=140 K
900 800 700
σ (a.u.)
600
T=360 K
1200 1000 800 600 0
20
40
60
80
t (ps)
Fig. 2. Time evolution of atomic equivalence indexes of K1Na54, averaged over 1000 time steps, for representative temperature values. The r curve corresponding to the K impurity is represented by a bold line.
177
are not useful to detect the second step in the melting transition. The reason is that K diffusion is restricted to the surface shell for all these temperatures. Segregation of K to the surface liquid layer is the expected behavior in terms of surface energy arguments. However, at T = 360 K, the impurity atom is able to make temporary excursions to the core atomic layers, as seen in Fig. 2, and mixing of the two atomic species is complete. Thus, some high-energy isomers (see Table 1) with an interior K atom are visited at this temperature, and this is reflected in a small peak in the specific heat. However, as the next higher temperature considered, evaporation of some atoms was observed, so that the small peak in the specific heat can also be a signature of incipient evaporation.
4. Summary In this paper, an efficient orbital-free molecular dynamics method has been employed to study the melting-like transition in K1Na55. The accuracy of the proposed scheme has been demonstrated by direct comparison to ab initio Kohn–Sham calculations. Our calculations predict icosahedral structures to be more stable than cuboctahedral structures. The lowest energy isomer is found with the K impurity located at a vertex surface position. The melting-like transition is shown to proceed stepwise, although the thermal evolution of the specific heat presents just one broad peak. Thus, we have analyzed a variety of melting indicators in order to obtain a complete picture of the melting mechanism. In a first stage, at T = 130 K, the impurity atom is the only one which does not participate in diffusion. This is probably due to the larger size of the K atom inhibiting diffusive motions inside a host cluster formed by Na atoms. In a second step, at T = 150 K, the K impurity begins to participate in diffusion events. K impurity is found to visit the cluster inner region only for sufficiently high temperatures (as high as 360 K, very close to the evaporation limit). Melting-like transitions in clusters are usually detected by examining the thermal evolution of the specific heat. This is a very important indicator as it can be experimentally determined [1a]. However, being a thermal (as opposed to structural) indicator, it does not provide enough information about the mechanisms by which melting proceeds. In this paper, we have shown the specific heat peak locates the beginning of the melting-like transition of K1Na54. However, at a temperature corresponding to the main peak in the specific heat, the cluster can not be considered to be fully liquid. The melting process in K1Na54 is peculiar enough that even the ‘‘atomic equivalence indexes’’, which provide very detailed structural information, are not able to
178
A. Aguado et al. / Computational Materials Science 35 (2006) 174–178
detect the second step in the transition. An analysis of the temperature dependence of diffusion constants for each atomic species, however, provides the dynamical information required to obtain a complete picture of the melting process in K1Na54. Doping of Na55 with K (and also with Li and Cs) [16] is found to significantly lower the melting temperature (as determined from the main peak in the specific heat). At the same time, however, due to incomplete mixing of different atomic species, the average distribution of atoms becomes homogeneous only at much higher temperatures, which shows the melting-like transition proceeds in a wider temperature range for these impurity-doped clusters. The impurity atom begins to participate in diffusion events at different temperatures than host (Na) atoms, which provides another mechanism for stepwise melting in multicomponent clusters, not present in one-component clusters. We will leave further investigations on melting of multicomponent clusters of different concentrations (nanoalloys) to future studies.
Acknowledgements This work was supported by Junta de Castilla y Leo´n (Project VA073/02) and DGES (Project MAT200204393-C02-01). A. Aguado also acknowledges financial support from the Spanish Ministry of Science and Technology, under the Ramo´n y Cajal program.
References [1] (a) M. Schmidt, R. Kusche, W. Kronmu¨ller, B. von Issendorff, H. Haberland, Phys. Rev. Lett. 79 (1997) 99; (b) M. Schmidt, R. Kusche, B. von Issendorff, H. Haberland, Nature 393 (1998) 238; (c) R. Kusche, Th. Hippler, M. Schmidt, B. von Issendorff, H. Haberland, Eur. Phys. J. D 9 (1999) 1; (d) M. Schmidt, H. Haberland, Compt. Rend. Physique 3 (2002) 327. [2] K. Manninen, A. Rytko¨nen, M. Manninen, Eur. Phys. J. D 29 (2004) 39. [3] (a) A. Aguado, J.M. Lo´pez, J.A. Alonso, M.J. Stott, J. Chem. Phys. 111 (1999) 6026; (b) A. Aguado, J.M. Lo´pez, J.A. Alonso, M.J. Stott, J. Phys. Chem. B 105 (2001) 2386. [4] Y.G. Chushak, L.S. Bartell, J. Phys. Chem. B 107 (2003) 3747. [5] M.J. Lo´pez, P.A. Marcos, J.A. Alonso, J. Chem. Phys. 104 (1996) 1056. [6] S. Huang, P.P. Balbuena, J. Phys. Chem. B 106 (2002) 7225. [7] K. Joshi, D.G. Kanhere, Phys. Rev. A 65 (2002) 043203. [8] K. Joshi, D.G. Kanhere, J. Chem. Phys. 119 (2003) 12301. [9] D.S. Mainardi, P.P. Balbuena, Int. J. Quantum Chem. 85 (2001) 580. [10] M. Pearson, E. Smargiassi, P.A. Madden, J. Phys. Condens. Mat. 5 (1993) 3221. [11] P. Hohenberg, W. Kohn, Phys. Rev. 136 (1964) 864B. [12] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) 1133A. [13] W.C. Swope, H.C. Andersen, J. Chem. Phys. 76 (1982) 637. [14] J.M. Soler, E. Artacho, J.D. Gale, A. Garcı´a, J. Junquera, P. Ordejo´n, D. Sa´nchez-Portal, J. Phys. Condens. Mat. 14 (2002) 2475. [15] V. Bonacic´-Koutecky´, J. Jellinek, M. Wiechert, P. Fantucci, J. Chem. Phys. 107 (1997) 6321. [16] A. Aguado, L.E. Gonza´lez, J.M. Lo´pez, J. Phys. Chem. B 108 (2004) 11722.
Molecular dynamics simulation of the melting-like transition in K1Na54 Andre´s Aguado
*,1,
Jose´ M. Lo´pez, Sara Nu´n˜ez
Department of Theoretical Physics, University of Valladolid, Valladolid 47011, Spain Received 24 April 2004; accepted 11 October 2004
Abstract A set of molecular dynamics simulations at constant total energy has been performed in order to investigate the melting-like transition of the impurity-doped cluster K1Na54. An orbital-free density-functional-theory technique is employed to obtain the forces on atoms in an efficient way. The total simulation time is approximately 1.5 ns, which is required to obtain statistically meaningful thermal averages. The presence of just one K impurity reduces considerably the melting temperature value (as compared to Na55). Although only one peak is observed in the thermal evolution of the specific heat, analysis of the diffusion constants shows that the cluster melts in two separate steps: Na atoms are able to diffuse at lower temperatures than the K impurity. 2005 Elsevier B.V. All rights reserved. PACS: 31.15.Qg; 36.40.Ei; 36.40.Sx; 64.70.Dv Keywords: Phase transitions in clusters; Alkali clusters; Ab initio molecular dynamics
1. Introduction A series of recent experiments by Haberland and coworkers [1a–d] has renewed the interest in the melting transition of small atomic clusters. These experiments reveal that variation of melting temperature Tm with cluster size N is non-monotonic, contrary to expectations based on classical arguments stating that melting temperature of a finite liquid drop should monotonically decrease as its radius shrinks. The detailed evolution of Tm with N shows enhanced stability of the solid phase for selected ranges of cluster sizes which are close to either electronic and/or geometrical shell closings. Thus, it is now accepted that both electronic and structural
*
Corresponding author. Tel.: +34 83 423147; fax: +34 83 423013. E-mail address: [email protected] (A. Aguado). 1 Supported by the ‘‘Ramo´n y Cajal’’ program (Spanish Ministry of Science and Technology). 0927-0256/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2004.10.011
effects will play a role in explaining the observed results, but apart from this tacit assumption, no definite explanation of the experimental observations has emerged yet. Computer simulation of the melting process may provide physical insight into the microscopic details of the melting-like transition and, hopefully, reproduce and explain the experimental Tm(N) curve. A very recent attempt in this direction has been given by Manninen et al. [2] These authors perform molecular dynamics (MD) simulations of melting in NaN (N = 40–355) and reproduce non-monotonic Tm(N) behavior. However, due to computational limitations, interatomic forces in all clusters with N P 93 are calculated employing a phenomenological potential, and the Tm(N) curve is not in good agreement with the experimental one for N P 142. The reasons for this discrepancy may be: (a) purely electronic effects (not well described by any phenomenological potential) have a strong influence; (b) interatomic forces provided by that potential might be not accurate
A. Aguado et al. / Computational Materials Science 35 (2006) 174–178
enough so that correct ground state (GS) isomers (for example, those predicted by an ab initio methodology) are not reproduced; (c) even if the potential correctly reproduces the ab initio energetic ordering of isomers, it is very difficult to locate the GS isomer when so many atoms are involved; etc. Agreement with experiment is much better for smaller clusters. Specifically, Manninen et al. reproduce the experimental result T m ðNaþ 55 Þ > T m ðNaþ 93 Þ, and MD simulations by Aguado et al. [3a,b] reproduce experimental melting temperatures of selected clusters in the range N = 8–142 to within a few percent. In the future, employment of more accurate interatomic forces and a thorough search for GS isomers may lead to satisfactory results also for large sizes. Consideration of the melting-like transition in impurity-doped metal clusters (and in multicomponent clusters in general) is even more complicate from a computational point of view. First of all, compositional ordering increases significantly the number of different isomer configurations, and thus the search for the GS isomer is much more time consuming. Also, surface segregation of the component with lowest surface energy may result in specific features such as partial melting of one component or incomplete mixing. Finally, and most important for the melting problem, the presence of impurities in an otherwise homogeneous cluster might alter structural properties, influencing thermal properties and, in particular, the location of the melting point. Although most of the work on multicomponent clusters has been devoted to structural properties, a few molecular dynamics simulations exist studying their freezing [4] and melting [5–8] transitions, and the influence of temperature on segregation [9]. In the present paper we report the results of MD simulations of the melting-like transition of the impuritydoped cluster K1Na54. To better identify the details of the solid-to-liquid transition, we have employed an orbital-free (OF) version [10,3a] of density-functional-theory (DFT) [11] (with core electrons replaced by pseudopotentials) in the calculation of interatomic forces, which are then used in a classical (Newtonian) description of atomic motions. Computational expense scales linearly with system size in OF approaches, allowing consideration of MD simulation runs long enough to extract meaningful statistical averages. The necessary employment of local pseudopotentials and approximate density functionals for the electronic kinetic energy reduces the accuracy of these methods, as compared to traditional Kohn–Sham (KS) approaches [12]. Nevertheless, OFMD includes an approximate explicit description of the electronic density, and thus provides a more transferable description of metal clusters than phenomenological potentials, for example. At the same time, a better statistical sampling of phase space (as compared to traditional KS approaches) may be achieved with only moderate effort. We will show that introduction of a
175
single K impurity in a Na55 host significantly affects thermal behavior: (1) the melting temperature is considerably lower than that of homogeneous Na55; (2) melting in K1Na54 proceeds in two separate steps, contrary to Na55, which melts in one single stage [3b]. We will also show the necessity of evaluating several melting indicators in order to obtain a complete picture of the melting process. In Section 2 we briefly present some technical details of the method. The results are presented and discussed in Section 3 and, finally, Section 4 summarizes our main conclusions.
2. Details of computer simulations Our implementation of OFMD is based on an expansion of the electron density in a basis set formed by plane waves. In order to make this implementation efficient, local pseudopotentials are employed to represent the effect of core electrons on the valence density. Thus, only this valence density is explicitly described. Details of the specific pseudopotentials and electronic kinetic and exchange-correlation energy functionals employed are the same as in our previous work [3a,3b]. The cluster is placed in a unit cell of a cubic superlattice, and the set of plane waves periodic in the superlattice is used as a basis set to expand the valence density. ˚ . and an Present calculations use a supercell of edge 35 A energy cut-off in the plane wave expansion of the density of 20 Ry. A 96 · 96 · 96 grid is used to perform Fast Fourier Transforms. The equations of motion are integrated using the velocity Verlet algorithm [13] with a time step equal to 4(3) fs for temperatures lower (higher) than 100 K. These choices resulted in a conservation of the total energy better than 0.1%. The ground state electron density is found for each configuration of atoms visited in the course of an MD trajectory, which allows the accurate calculation of atomic forces by means of the Hellman–Feynman theorem. As a first step, we have tried to locate a plausible GS isomer for K1Na54. Specifically, we have considered complete icosahedral and cuboctahedral structures and all inequivalent positions for the substitutional impurity. In order to explicitly check the reliability of OF predictions, we optimize the same structures employing pseudopotential KS-DFT calculations as implemented in SIESTA [14], applying the same exchange-correlation functional (LDA) for consistency. Standard double-zeta plus polarization (DZP) basis sets were adopted for each atom, together with an energy cutoff of 20 Ry. Stable isomer geometries are found by employing a conjugate gradients optimization technique and requiring the residual force on each atom to be smaller than ˚. 0.001 eV/A Once we have found a candidate GS isomer, we perform a set of constant energy MD simulations in order
A. Aguado et al. / Computational Materials Science 35 (2006) 174–178
Different values of the average temperature were obtained by isokinetic thermalization runs (approximately 5 ps long) previous to each constant energy run. A typical simulation time for the constant energy runs of solid-like clusters is 50 ps, but runs longer than 200 ps were performed for those energies where the melting transition sets in. For energies corresponding to liquidlike clusters, the length of the simulations was always in the range 70–100 ps. The total simulated time for each cluster was approximately 1.5 ns.
3. Results There are four geometrically inequivalent sites in a perfect 55-atom icosahedron. Introducing a substitutional impurity in each of these sites generates four different K1Na54 isomers, whose energies are compared in Table 1. We do not include results for corresponding cuboctahedral isomers as these are located at higher energies. As expected from size mismatch and surface energy arguments, K prefers to be located at a vertex site of the cluster surface. Moreover, OF calculations reproduce the energetic ordering of isomers predicted by SIESTA, which supports the employment of our OF procedure to study thermal properties. From a structural point of view, the only appreciable difference between OF and SIESTA calculations is that OF distances are systematically longer by slightly less than 1%, an error which lies well within the expected accuracy of our methodology. In Fig. 1, we show the caloric curve, specific heat, partial rms bond length fluctuations and diffusion constants of K1Na54, as a function of average temperature. A thermal phase transition is indicated in the caloric
Table 1 Energies (in meV/atom) of several isomers of K1Na54, taking the energy of the corresponding ground state (GS) isomer as a reference
SIESTA OF
a
b
c
d
9.5 10.0
7.6 8.0
GS GS
1.5 1.4
K atom is located at (a) the center of the icosahedron; (b) any of the 12 equivalent positions within the first icosahedral shell; (c) a vertex position of the surface shell (again, the multiplicity of this site is 12); (d) any of the 30 non-vertex sites of the outer shell.
-12.4 -12.45 Na-Na Na-K
-12.5 -12.55 7
0.2 0.15 0.1
δ
0.3 0.25
0.05 0 15
Na K
6
10
5 4
5
D (bohr2/ps)
j
-12.3 -12.35
E (hartree)
to obtain the caloric curve, specific heat, root-meansquared bond length fluctuations and diffusion constants. The way these averages are evaluated is explained in previous works [3a,b]. We just show here explicitly the expression for the ‘‘atomic equivalence indexes’’ [15], as these are less standard quantities X ri ðtÞ ¼ j~ Ri ðtÞ ~ Rj ðtÞj. ð1Þ
Specific Heat (kB)
176
3 0
2 0
100
200
T (K)
300
400
100
200
300
400
T (K)
Fig. 1. Caloric and specific heat curves (left side) and partial rms bond length fluctuations and diffusion constants (right side) of K1Na54, taking the internal cluster temperature as the independent variable.
curve by a change of slope, the slope being the specific heat; the height of the step gives an estimate of the latent heat of fusion. However, melting processes are more easily recognized as peaks in the specific heat or steps in d, as a function of temperature. All three indicators lead to the same value for the estimated melting temperature (130 K), which is a good check for the convergence of our results. Our estimation for the latent heat of fusion is of 0.35 eV. The temperature evolution of the specific heat displays just one main peak centered at Tm = 130 K. At this temperature, there is an abrupt increase in the phase space volume available to the system. No more features are observed in the specific heat, except for temperatures close to 400 K, where evaporation of atoms sets in (incidentally, this last temperature is quite close to the evaporation temperature estimated by Schmidt and Haberland [1d] for Naþ 139 . An examination of the boiling transition is, however, outside the scope of the present communication). Analysis of the temperature evolution of diffusion constants, however, shows that the K atom does not participate in diffusion events until a temperature of 150 K is reached. This is significantly higher than Tm and shows that melting proceeds in two well defined steps: (1) at T = Tm, all sodium atoms (both surface and inner-shell atoms) can freely diffuse across the whole cluster volume, but K atom stays in its vertex position at the surface, without participating in diffusion processes; (2) at T = 150 K, free energy barriers against K diffusion may be overcome and the whole cluster is fluid-like. This two-step melting process is validated by visual inspection of the movement of atoms on a computer screen. The second step is not reflected in the specific heat because the relative increase in the available phase space associated to it is very small. It is not reflected either in the partial bond length fluctuations because diffusion of only Na atoms is enough to induce a step in the thermal evolution of Na–K bond length
A. Aguado et al. / Computational Materials Science 35 (2006) 174–178
fluctuations. It only can be appreciated in indicators which depend on the identity of each single atom, as is the case for the diffusion constants. Fig. 2 shows the time evolution of short-time-averages of the atomic equivalence indexes introduced in previous section, for several representative temperatures. These quantities are adequate to analyze some microscopic details of the melting transition. At T = 120 K, atoms are grouped in shells, representing symmetry equivalent atomic sites. For the case at hand, we find the smallest value of r for the Na atom at the center of the icosahedron; then a group of 12 atoms forming the first icosahedral shell; next, a group of 30 lines representing those atoms of the outer shell capping the edges of the inner Na13 icosahedron; and finally, 12 lines for the atoms at vertex positions of the surface. Within this surface shell, K is clearly identified due to its larger nearest neighbor distances. Although no mixing between different shells is observed, some intrashell disorder already appears at this temperature which is very close to Tm. This disorder is induced by temporary excursions of some atoms to positions of low average radial atomic density. After these temporary excursions, each atom comes back to its original site, that is, this disorder does not involve intrashell diffusion. At a higher temperature (T = 140 K), the amplitude of these motions is large enough to induce mixing between different shells and melting of the whole ‘‘Na sublattice’’. Note that there is not a separate surface-melting stage, as observed for example in Na20 or Na92 [3a,b], but that both surface and inner Na atoms start to diffuse at the same temperature. This is similar to the homogeneous melting process observed for Na55 [3b]. The only difference is that, at this temperature, K atom is not able to diffuse yet. For all temperatures between 140 and 360 K, the atomic equivalence indexes look quite similar, that is, these indicators
σ (a.u.)
1000 900 800
T=120 K
σ (a.u.)
700 600 1000
T=140 K
900 800 700
σ (a.u.)
600
T=360 K
1200 1000 800 600 0
20
40
60
80
t (ps)
Fig. 2. Time evolution of atomic equivalence indexes of K1Na54, averaged over 1000 time steps, for representative temperature values. The r curve corresponding to the K impurity is represented by a bold line.
177
are not useful to detect the second step in the melting transition. The reason is that K diffusion is restricted to the surface shell for all these temperatures. Segregation of K to the surface liquid layer is the expected behavior in terms of surface energy arguments. However, at T = 360 K, the impurity atom is able to make temporary excursions to the core atomic layers, as seen in Fig. 2, and mixing of the two atomic species is complete. Thus, some high-energy isomers (see Table 1) with an interior K atom are visited at this temperature, and this is reflected in a small peak in the specific heat. However, as the next higher temperature considered, evaporation of some atoms was observed, so that the small peak in the specific heat can also be a signature of incipient evaporation.
4. Summary In this paper, an efficient orbital-free molecular dynamics method has been employed to study the melting-like transition in K1Na55. The accuracy of the proposed scheme has been demonstrated by direct comparison to ab initio Kohn–Sham calculations. Our calculations predict icosahedral structures to be more stable than cuboctahedral structures. The lowest energy isomer is found with the K impurity located at a vertex surface position. The melting-like transition is shown to proceed stepwise, although the thermal evolution of the specific heat presents just one broad peak. Thus, we have analyzed a variety of melting indicators in order to obtain a complete picture of the melting mechanism. In a first stage, at T = 130 K, the impurity atom is the only one which does not participate in diffusion. This is probably due to the larger size of the K atom inhibiting diffusive motions inside a host cluster formed by Na atoms. In a second step, at T = 150 K, the K impurity begins to participate in diffusion events. K impurity is found to visit the cluster inner region only for sufficiently high temperatures (as high as 360 K, very close to the evaporation limit). Melting-like transitions in clusters are usually detected by examining the thermal evolution of the specific heat. This is a very important indicator as it can be experimentally determined [1a]. However, being a thermal (as opposed to structural) indicator, it does not provide enough information about the mechanisms by which melting proceeds. In this paper, we have shown the specific heat peak locates the beginning of the melting-like transition of K1Na54. However, at a temperature corresponding to the main peak in the specific heat, the cluster can not be considered to be fully liquid. The melting process in K1Na54 is peculiar enough that even the ‘‘atomic equivalence indexes’’, which provide very detailed structural information, are not able to
178
A. Aguado et al. / Computational Materials Science 35 (2006) 174–178
detect the second step in the transition. An analysis of the temperature dependence of diffusion constants for each atomic species, however, provides the dynamical information required to obtain a complete picture of the melting process in K1Na54. Doping of Na55 with K (and also with Li and Cs) [16] is found to significantly lower the melting temperature (as determined from the main peak in the specific heat). At the same time, however, due to incomplete mixing of different atomic species, the average distribution of atoms becomes homogeneous only at much higher temperatures, which shows the melting-like transition proceeds in a wider temperature range for these impurity-doped clusters. The impurity atom begins to participate in diffusion events at different temperatures than host (Na) atoms, which provides another mechanism for stepwise melting in multicomponent clusters, not present in one-component clusters. We will leave further investigations on melting of multicomponent clusters of different concentrations (nanoalloys) to future studies.
Acknowledgements This work was supported by Junta de Castilla y Leo´n (Project VA073/02) and DGES (Project MAT200204393-C02-01). A. Aguado also acknowledges financial support from the Spanish Ministry of Science and Technology, under the Ramo´n y Cajal program.
References [1] (a) M. Schmidt, R. Kusche, W. Kronmu¨ller, B. von Issendorff, H. Haberland, Phys. Rev. Lett. 79 (1997) 99; (b) M. Schmidt, R. Kusche, B. von Issendorff, H. Haberland, Nature 393 (1998) 238; (c) R. Kusche, Th. Hippler, M. Schmidt, B. von Issendorff, H. Haberland, Eur. Phys. J. D 9 (1999) 1; (d) M. Schmidt, H. Haberland, Compt. Rend. Physique 3 (2002) 327. [2] K. Manninen, A. Rytko¨nen, M. Manninen, Eur. Phys. J. D 29 (2004) 39. [3] (a) A. Aguado, J.M. Lo´pez, J.A. Alonso, M.J. Stott, J. Chem. Phys. 111 (1999) 6026; (b) A. Aguado, J.M. Lo´pez, J.A. Alonso, M.J. Stott, J. Phys. Chem. B 105 (2001) 2386. [4] Y.G. Chushak, L.S. Bartell, J. Phys. Chem. B 107 (2003) 3747. [5] M.J. Lo´pez, P.A. Marcos, J.A. Alonso, J. Chem. Phys. 104 (1996) 1056. [6] S. Huang, P.P. Balbuena, J. Phys. Chem. B 106 (2002) 7225. [7] K. Joshi, D.G. Kanhere, Phys. Rev. A 65 (2002) 043203. [8] K. Joshi, D.G. Kanhere, J. Chem. Phys. 119 (2003) 12301. [9] D.S. Mainardi, P.P. Balbuena, Int. J. Quantum Chem. 85 (2001) 580. [10] M. Pearson, E. Smargiassi, P.A. Madden, J. Phys. Condens. Mat. 5 (1993) 3221. [11] P. Hohenberg, W. Kohn, Phys. Rev. 136 (1964) 864B. [12] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) 1133A. [13] W.C. Swope, H.C. Andersen, J. Chem. Phys. 76 (1982) 637. [14] J.M. Soler, E. Artacho, J.D. Gale, A. Garcı´a, J. Junquera, P. Ordejo´n, D. Sa´nchez-Portal, J. Phys. Condens. Mat. 14 (2002) 2475. [15] V. Bonacic´-Koutecky´, J. Jellinek, M. Wiechert, P. Fantucci, J. Chem. Phys. 107 (1997) 6321. [16] A. Aguado, L.E. Gonza´lez, J.M. Lo´pez, J. Phys. Chem. B 108 (2004) 11722.