Evaporation of a mixed cluster: KrXe13

Chemical Physics Letters 381 (2003) 471–478 www.elsevier.com/locate/cplett

Evaporation of a mixed cluster: KrXe13 P. Parneix b

a,*

, Ph. Br
echignac a, F. Calvo

b

a Laboratoire de Photophysique Mol
eculaire, B^ at. 210, Universit
e Paris-Sud, F91405 Orsay Cedex, France Laboratoire de Physique Quantique, IRSAMC, Universit
e Paul Sabatier, 118 Route de Narbonne, F31062 Toulouse Cedex, France

Received 22 September 2003; in final form 22 September 2003 Published online: 24 October 2003

Abstract The competition between two cluster fragments (KrXe12 and Xe13 ) in the evaporative process of the non-rotating mixed Lennard–Jones cluster KrXe13 has been investigated both from molecular dynamics simulations and Phase Space Theory (PST). The branching ratio between the two fragments is seen to be strongly altered by the thermodynamical changes, especially the preliminary core/surface isomerization of the Kr atom. Such a multistep melting behavior is well reproduced by PST when numerically exact radial potentials along the dissociative coordinate and anharmonic vibrational densities of states are considered. Ó 2003 Elsevier B.V. All rights reserved.

1. Introduction The mechanisms of cluster evaporation have attracted a significant attention from the experimental and theoretical communities [1–13]. While most of these works have been concerned with homogeneous systems, mixed clusters are intrinsically good candidates for studying the competition between different channels in the evaporation process, because of the various possible dissociation channels. Such an analysis is of particular importance in cluster physics. As an example, it is wellknown that a competition between monomer and dimer occurs in metallic clusters [1,2] but also in

*

Corresponding author. Fax: +33-169156777. E-mail address: [email protected] (P. Parneix).

small carbon clusters [3] as well as charged rare gas clusters [4]. Experimental information about these dissociation events is always obtained in a microsecond time domain, which corresponds to a relatively low temperature of the parent cluster, as compared to what can be achieved through bruteforce molecular dynamics computer simulations. It is thus especially important to develop a theoretical approach alternative to simulation. These last years many researchers have used different statistical theories to describe evaporation in clusters [5–13]. From a proper comparison with benchmark MD simulations, Phase Space Theory (PST) [14–16] now clearly appears as the most quantitatively appropriate for describing cluster evaporation. A particularly interesting feature shown by these theoretical studies is the sensitivity of the

0009-2614/$ - see front matter Ó 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2003.10.005

472

P. Parneix et al. / Chemical Physics Letters 381 (2003) 471–478

mean kinetic energy release (KER) with respect to the thermodynamical state of the product [6,8]. One could thus imagine using evaporation dynamics as an indirect probe of the solid–liquid phase change in clusters, also well characterized recently in calorimetry type experiments [17]. Most of these theoretical studies focused on the evaluation of the evaporation rate and the KER distribution in a non-rotating cluster. However, it was shown that rotating systems could be treated as well, [11–13] which should not be surprising since conservation of angular momentum J is rigorously included in PST. In the present work, we mainly focus on the calculation of the branching ratio in the evaporation of a non-rotating mixed cluster, namely KrXe13 . Since we are dealing with two possible products, the connection between the observed KER and the thermodynamical characteristics may not be as straightforward as in the case of an homogeneous system. The hope is to be able to use the branching ratio as an extra indicator of the possible isomerization processes. Recently, [18] we have followed the evolution of the branching ratio in a mixed Ar4 X3 cluster described using Lennard–Jones (LJ) interactions. X was there an arbitrary model atom, characterized by its reduced LJ parameters of distance, energy and mass r
¼ rX–X =rAr–Ar , e
¼ eX–X =eAr–Ar , and M
¼MX =MAr . It was noticed that the smaller r
 M
 e
, the larger the probability for evaporating a X-type atom rather than an Ar atom [18]. If r
M
e
is too small, the branching ratio is always equal to 0 or 1. On the other hand, if r
M
e
is close to 1, the dynamics becomes slower and it is very hard to simulate a large set of evaporative trajectories close to the melting energy. Considering these practical constraints suggests the best choice for realistic atomic LJ clusters to be the Kr/Xe system for which r
M
e
¼ 0:434. More precisely, the dynamics of the parent KrXe13 cluster has been investigated in details. In the next section, we briefly recall the main equations underlying the PST calculation of the mean KER and of the probability to have a specific fragment after evaporation. In Section 3 we discuss our results before finally concluding.

2. Methodology We consider a parent cluster Krn Xem , also denoted simply as ðn; mÞ in what follows, characterized by a total vibrational energy E and a total angular momentum J ¼ 0. After evaporation of a monomer, two product clusters can be formed, namely Krn1 Xem or Krn Xem1 . Following Weerasinghe and Amar [6] and Jarrold [19], the dissociation statistics in a given channel characterized by its couple of integers ðn0 ; m0 Þ is governed by the differential rate given by: Rn0 ;m0 ðetr ; E; J ¼ 0Þ ðn0 ;m0 Þ

¼ R0

Xn0 ;m0 ðE  E0  etr ÞCn0 ;m0 ðetr Þ : Xn;m ðEÞ

ð1Þ

In the latter equation, R0 is a constant factor that accounts for channel and rotational degeneracies [6], Xn;m and Xn0 ;m0 are the vibrational densities of states (VDOS) of the parent Krn Xem and product Krn0 Xem0 clusters, respectively. Cn0 ;m0 is the rotational density of states (RDOS) of the fragments. ðn0 ;m0 Þ Finally E0 is the energy difference between the most stable configurations of the parent and product clusters, and etr corresponds to the total kinetic energy released in the evaporation process. Obviously two differential rates have to be considered for a mixed cluster. The first one corresponds to the loss of a Kr atom ðn0 ¼ n  1; m0 ¼ mÞ, the second being the loss of a Xe atom ðn0 ¼ n; m0 ¼ m  1Þ. From the differential rate, the probability for finding etr between etr and etr þ detr in a given channel is: Pn0 ;m0 ðetr ; E; J ¼ 0Þ ¼ Rn0 ;m0 ðetr ; E; J ¼ 0Þ
Z

ðn0 ;m0 Þ

EE0

Rn0 ;m0 ðetr ; E; J : ¼ 0Þdetr :

ð2Þ

0

The knowledge of the differential rate function Rn0 ;m0 readily leads to the average kinetic energy released: hetr i

ðn0 ;m0 Þ

¼

Z 0

ðn0 ;m0 Þ

EE0

etr Pn0 ;m0 ðetr ; E; J ¼ 0Þdetr :

ð3Þ

P. Parneix et al. / Chemical Physics Letters 381 (2003) 471–478

In the PST framework, the ratio between the probability to evaporate a Xe atom to that for a Kr atom is proportional to the ratio between two integrals involving the VDOS of these two product clusters. Following previous works [8,18], this ratio, denoted as c, can be expressed as: R EE0ðn1;mÞ ðn1;mÞ Cn1;m ðetr ÞXn1;m ðE  E0  etr Þdetr c ¼ R0 ðn;m1Þ : EE0 ðn;m1Þ C ðe ÞX ðE  E  e Þde n;m1 tr n;m1 tr tr 0 0 ð4Þ Finally the percentage of evaporated Kr atom, denoted as PKr , is directly given by: c : ð5Þ PKr ¼ 100  1þc ðn0 ;m0 Þ

The calculation of hetr i and PKr requires the computation of Xn0 ;m0 and Cn0 ;m0 for both product/ clusters, but the VDOS of the parent cluster is not needed. In the following, we briefly indicate how the rotational and vibrational DOS were calculated. First, we focus on the RDOS. By taking into account the mechanical constraints due to the conservation laws (energy and angular momentum) in the dissociation, Chesnavich and Bowers 2 [16] have shown that Cn0 ;m0 is equal to Jmax if the product cluster is considered as a spherical rigid body. Here Jmax is the largest possible value of the angular momentum of the product cluster. Then C is obtained from an extra relationship between etr 2 and Jmax , which follows from the expression [6]. 0

0

2 ¼ etr ; ey ðJmax Þ þ Bsðn ;m Þ Jmax

ð6Þ

ey being the centrifugal barrier for the orbital angular momentum L ¼ Jmax . In the case of a radial dissociation potential with the form V ðrÞ ¼ C6 =r6 , ey takes the simple

473

3 =Kn0 ;m0 with Kn0 ;m0 given by analytical form Jmax 3=2

1=2

Kn0 ;m0 ¼ ð6ln0 ;m0 Þ C6 =2. ln0 ;m0 , the reduced mass of the evaporating cluster, is reported in Table 1 0 0 for the two exit channels. Bsðn ;m Þ corresponds to the rotational constant of the product cluster Krn0 Xem0 considered here as a spherical top. As shown previously [6,11], the radial potential is poorly reproduced by a simple C6 =r6 law. The exact radial potential between the atomic fragment and the product cluster has been calculated using the Wang–Landau method [20], adapted for calculating potential of mean forces and average potentials [21]. The interactions between the fragments have thus been fitted to C6 =ðr  r0 Þ6 laws, with parameters given in Table 1. From these radial potentials, ey can be easily deduced as a function of the orbital angular momentum. As concerns the VDOS, distributions of potential energy at various temperatures have been extracted from parallel tempering Monte Carlo simulations [22]. From these distributions, the configurational DOS was built from the multiple histogram method [23] and the VDOS was obtained by a convolution with the analytical kinetic energy DOS. Finally the evaporation dynamics has been explored from MD simulations using a set of 10 000 independent trajectories. A fifth-order Adams–Moulton predictor– corrector algorithm has been used to propagate evaporative trajectories from which the mean kinetic energies released and the branching ratio were computed. More computational details have been given in a previous article [8]. As the number of Xe atoms is significantly larger than Kr atoms, note that all the results will be given in LJ units corresponding to Xenon,  for distances, m ¼ 131:29 a.m.u. that is r ¼ 4:1 A for masses, and e ¼ 221 K for energies and temperatures.

Table 1 Binding energies, interaction parameters, rotational constants and reduced masses of the product clusters and corresponding fragments from the evaporation of KrXe13 Product cluster

Emin ðeÞ

C6 ðer6 Þ

r0 ðrÞ

BT ¼0 ðmr2 Þ

lðmÞ

Xe13 KrXe12

)44.33 )44.15

20.372 33.614

0.506 0.493

0.0534 0.0559

0.6084 0.9268

All quantities are expressed in LJ units of the Xe atom. At the (common) melting point, the rotational constants have been scaled by 0.85 from the 0 K value.

474

P. Parneix et al. / Chemical Physics Letters 381 (2003) 471–478

3. Results and discussion Before discussing the results of the evaporation dynamics, we give the main features of the two product clusters involved in the dissociation. Not surprisingly, the minimum energy structures of the two product clusters are icosahedral. In the case of KrXe12 , the Kr atom is located at the center of the cluster. This may not have been obvious, because Kr–Xe bonds are weaker than Xe–Xe bonds, in addition to being shorter. Actually this structure results from the compression of the bonds between the central and surface atoms, which favors shorter bonds, hence a Krypton core. The values of the binding energies associated to the global minima of KrXe12 and Xe13 are reported in Table 1. Their difference is only 0.18. Since we are interested in the sensitivity of the branching ratio to the solid–liquid phase change, we have computed some thermodynamical properties of these two clusters. In Fig. 1 the heat capacities have been plotted as a function of temperature. In the mixed KrXe12 cluster, melting as indicated by the top of the peak appears at a slightly smaller temperature than in the neat Xenon cluster. An another interesting feature in the caloric curve is the presence of a smaller peak around T ¼ 0:14 for KrXe12 . To interpret this feature, we have com2 puted the average square distance hdKr i between the Kr atom and the cluster center-of-mass. The variations of this quantity with temperature are plotted in Fig. 2. They show two successive in-

Heat capacity (kB)

KrXe12 Xe13

Fig. 2. Variations of the thermally average square distance 2 hdKr i between the Kr atom and the center of mass of KrXe12 as a function of temperature. Also superimposed is the probability P1 to find the global minimum, calculated in the harmonic superposition approximation with the two lowest isomers only. Both scales are in LJ units.

creases, at T ¼ 0:14 and 0.27, respectively, each of them being correlated with a bump in the heat capacity. The feature around T ¼ 0:14 therefore indicates that the Kr atom can jump from the center to the surface of the cluster, before even leaving the icosahedral shell at the melting point. The interpretation of the preliminary feature in the caloric curve as a core/surface isomerization has been further addressed using a simple twostates, harmonic superposition approximation [24]. In this model, the thermodynamics is built only from the two lowest isomers of the KrXe12 system, which turn out to be the two non-equivalent icosahedral isomers where the Kr atom occupies either the core site (Ih symmetry) or one surface site (C5v symmetry). The knowledge of the potential energies Ei and the geometrical averages  i of the vibrational frequencies of these isomers x allows to estimate the canonical probabilities Pi ðT Þ of finding each of them as the ratio of partition functions, Pi ðT Þ ¼ Qi ðT Þ=ðQ1 ðT Þ þ Q2 ðT ÞÞ with Qi ðT Þ ¼ ni

20

Temperature

Fig. 1. Heat capacities of KrXe12 and Xe13 .

ebEi :  i Þj ðbhx

ð7Þ

In the above equation, ni accounts for all permutational isomers having the same contribution and is given for the KrXeN system by ni ¼ 2N !=hi , with hi the order of the point group of structure i. b ¼ 1=kB T is the usual inverse temperature, and

P. Parneix et al. / Chemical Physics Letters 381 (2003) 471–478

j ¼ 3N  3 is the number of independent degrees of freedom. The variations of the probability P1 ðT Þ with temperature have been superimposed on the graph of Fig. 2. Even in the harmonic approximation, the sharp drop of P1 ðT Þ close to T ¼ 0:15 is fully 2 consistent with the observed increased in hdKr i. As explained in the previous section, the PST model has been used in the sphere þ atom approximation. One point to check now is the influence of the impurity (Kr atom) on the sphericity of the cluster as a function of temperature. To obtain this information, we have computed the rotational constants A, B and C from the Monte Carlo simulations. We note f the asphericity factor, defined as f¼

AC AC ; ¼3 Bs AþBþC

ð8Þ

where Bs ¼ ðA þ B þ CÞ=3 is the rotational constant of the associated spherical top. In Fig. 3(a), f has been plotted for the two clusters. f strongly increases for the two systems near melting, up to about 0.25 at T ¼ 0:28. However, this value reached by f remains relatively low with respect to very assymetric structures, for which f would have a magnitude closer to 1 or

475

even larger. Thus the spherical approximation seems reasonable, and the influence of the Kr atom appears relatively small. As a side note, the core/ surface isomerization exhibited by the Kr atom does not significantly affect asphericity, only a small increase in f being observed near T ¼ 0:14. In Fig. 3(b), the thermally averaged value of Bs has been plotted as a function of temperature. The values at T ¼ 0 are also reported in Table 1. These values are needed to calculate the RDOS. Since we are interested in the energy range for which the products are formed near the melting point, we will use values of the rotational constant calculated at T  0:28. The percentage of evaporated Kr atom, PKr , is plotted in Fig. 4 as a function of the total energy, from MD and PST for comparison. Because the anharmonic VDOS are obtained up to some multiplicative factor, we cannot directly obtain the exact value of PKr . Thus we have scaled PKr to reproduce the MD values at high energies. Within this approximation, it has to be noted that the PST calculation perfectly reproduces the energy dependence of PKr in a broad energy range. The PST formalism has also been used within the harmonic (or Engelking) approximation. As the harmonic VDOS are totally determined by the harmonic frequencies of the global minimum only, no

(a) 0.4

KrXe12 Xe13

100

ζ

0.3 0.2

80

0

0

0.1

0.2

0.3

0.4

PKr

0.1

60

(b) 0.06 40

Bs

0.05

KrXe12 Xe13

0.04 0.03

0

0.1

20 -45 0.2 Temperature

0.3

PST; harmonic PST; HSA PST; anharmonic MD -40

-35 Energy

-30

-25

0.4

Fig. 3. (a) Average asphericity factor f versus temperature for KrXe12 and Xe13 and (b) average rotational constant Bs versus temperature for the same clusters.

Fig. 4. Relative probability of evaporating a Krypton atom versus total energy (in LJ units) in several approximations. HSA refers to the harmonic superposition approximation to the densities of states based on the two lowest isomers of KrXe12 only.

476

P. Parneix et al. / Chemical Physics Letters 381 (2003) 471–478

further adjustment is required. The results of this basic harmonic approximation, also plotted in Fig. 4, badly fail in reproducing MD data. The disagreement between the harmonic and anharmonic descriptions starts clearly at around E ¼ 39, which corresponds to the onset of core/ surface isomerization of the Kr atom. Therefore it seems important to account for anharmonicities in the statistical treatment of cluster evaporation. This result confirms previous observations by Weerasinghe and amar [6], by Peslherbe and Hase [7,9] and more recently by us [8,11]. As can be noted in Fig. 4, a strong decrease in PKr is observed around this energy due to the large increase in the VDOS of KrXe12 near this preliminary isomerization. Since only the denominator in Eq. (4) is altered near the phase change, it explains why the effect on PKr is so strong. To emphasize the role of the second isomer of KrXe12 in the evaporation statistics, we have also calculated the branching ratio by assuming that only the two icosahedral isomers contribute to the VDOS. This is again a harmonic superposition approximation, which reads [24]: XðEÞ ¼ n1

ðE  E1 Þj1 ðE  E2 Þj1 þ n ; 2  1 Þj  2 Þj ð hx ð hx

RDOS is approximated as proportional to etr also gives a probability of emitting the Krypton atom in good agreement with MD. This situation corresponds to the neglect of the centrifugal barrier. However, it must be kept in mind that the Weisskopf theory is not expected to work straightforwardly in the more complicated cases of the dissociation of a dimer (or a larger fragment), because the RDOS no longer grows linearly with etr . More realistic treatments of PST, as in the sphere þ linear approximation, should then be considered instead [13]. Finally, we have represented in Fig. 5 the variations of the average kinetic energy released in the two dissociations processes KrXe13 ! Xe13 þ Kr (upper panel) and KrXe13 ! KrXe12 þ Xe (lower panel). This quantity is also known to be a sensitive probe of possible phase changes in the product cluster [6], which is also clearly apparent here for both channels of evaporation. As the total energy increases to E  32, the average KER reaches a plateau, indicative of melting in both products. The variations of hetr iðEÞ are made even clearer in

ð9Þ

provided that each of both terms in the right hand side are defined and positive. The branching ratio using this approximation for the VDOS of KrXe12 is also plotted in Fig. 4. Now the observed variations are much closer to the MD predictions, and also to the fully anharmonic PST calculation. This better agreement simply comes from the much larger statistical weight of the isomer with Kr at the surface, as previously seen in Fig. 2. Hence the core/surface preliminary isomerization plays here a very strong role in determining the products after evaporation. This effect is further reinforced by the vicinity of the melting points in the two clusters Xe13 and KrXe12 , because the two densities of vibrational states become then more similar at moderate energies E > 35. In addition, the influence of the rotational constant on the PST results has been considered. Taking for Bs the value at T ¼ 0 yields only a small difference with the value at the melting point. The Weisskopf-type [25] approximation in which the

Fig. 5. Average kinetic energy released upon evaporation of (a) a Kr atom (a) or (b) a Xe atom from Kr13 , as a function of total energy E. The dashed straight lines correspond to PST estimates within the simple harmonic approximation (Engelking approximation). The vertical arrow indicates an inflection due to core/surface isomerization. In each panel, the inset graph shows the derivative of hetr i with respect to E.

P. Parneix et al. / Chemical Physics Letters 381 (2003) 471–478

the plot of their derivative with respect to total energy, represented in the inset graphs. Nevertheless, an interesting difference can be seen between the two curves in Fig. 5. When E gets close to about )37 LJ units, there is a small but visible inflection in hetr i for the loss of a Xe atom, better seen in the bump in the derivative. This is precisely the energy at which the core/surface isomerization takes place. Therefore the average KER seems to contain as much information as the more usual thermodynamical indicators in providing evidences for phase changes in atomic clusters, and we find here some evidence that the variations of hetr i with energy closely match the shape of the caloric curve. To our knowledge, this is the first time that a multistep melting process is seen to be fully characterized by the evaporation statistics.

4. Conclusion The possibility to perfectly reproduce from a PST calculation the branching ratio between two different evaporative channels implying loss of monomers is very encouraging. The next natural step would be to directly describe the competition between monomer and dimer fragments using PST. For the dimer evaporation channel, the PST formalism recently developed for molecular systems should be used [13]. This could be achieved on Naþ n clusters, for which experimental data have been collected [1]. Alternatively, acetonitrile molecular clusters could provide a further candidate for PST calculations. These clusters have already been simulated by one of us [26], and their statistical evaporation was shown to display a monomer/dimer competition, apparently correlated with the phase change in the product cluster. Another interesting point of the present work is that we have clearly shown that a multistep melting process could be fully characterized by the evaporation statistics. Not only the branching ratio and the average KER are sensitive to the solid–liquid phase transition, but the core/surface preliminary isomerization also induces visible changes on these observables. It is interesting to

477

note that the branching ratio is a particularly well-suited observable to characterize this ‘‘premelting’’ dynamical behavior. Due to high barriers on the potential energy surface, this isomerization process will likely be very slow at low energies, and it should significantly hinder the evaporation of the Kr atom. Therefore the absolute rate constant for the dissociation of a Xenon atom should also be expected to undergo a drastic reduction as the total energy becomes lower than the isomerization barrier between the two lowest structures of KrXe12 . Such predictions could not be checked using brute force molecular dynamics, which is precisely why we should validate and use suitable statistical approaches such as PST.

Acknowledgements The authors wish to thank the GDR Agr
egats, Dynamique et R
eactivit
e for financial support.

References [1] C. Br
echignac, Ph. Cahuzac, J. Leygnier, J. Weiner, J. Chem. Phys. 90 (1989) 1492. [2] M. Vogel, K. Hansen, A. Herlert, L. Schweikhard, Appl. Phys. B 73 (2001) 411. [3] G. Martinet, M. Chabot, K. Wohrer, S. Della Negra, D. Gardes, J.A. Scarpaci, P. D
esesquelles, V. Lima, S. D
ıazTendero, M. Alcam
ı, P.-A. Hervieux, M.F. Politis, J. Hanssen, F. Mart
ın, Eur. Phys. J. D 24 (2003) 149. [4] see F. Calvo, J. Galindez, F.-X. Gad
ea, Phys. Chem. Chem. Phys. 5 (2003) 321, and references therein. [5] M.J. L
opez, J. Jellinek, Phys. Rev. A 50 (1994) 1445. [6] S. Weerasinghe, F.G. Amar, Z. Phys. D: At., Mol. Clusters 20 (1991) 167; J. Chem. Phys. 98 (1993) 4967. [7] G.H. Peslherbe, W.L. Hase, J. Chem. Phys. 105 (1996) 7432. [8] P. Parneix, F.G. Amar, Ph. Br
echignac, J. Phys. Chem. 239 (1998) 121. [9] G.H. Peslherbe, W.L. Hase, J. Phys. Chem. A 104 (2000) 10556. [10] P.A. Hervieux, B. Zarour, J. Hanssen, M.F. Politis, F. Mart
ın, J. Phys. B 34 (2001) 3331. [11] F. Calvo, P. Parneix, J. Chem. Phys. 119 (2003) 256. [12] P. Parneix, F. Calvo, J. Chem. Phys. 119 (2003) 0000. [13] F. Calvo, P. Parneix, J. Chem. Phys., submitted for publication. [14] P. Pechukas, J.C. Light, J. Chem. Phys. 42 (1965) 3281. [15] C.E. Klots, J. Phys. Chem. 75 (1971) 1526.

478

P. Parneix et al. / Chemical Physics Letters 381 (2003) 471–478

[16] W.J. Chesnavich, M.T. Bowers, J. Chem. Phys. 66 (1977) 2306. [17] M. Schmidt, T. Hippler, J. Donger, W. Kronm€ uller, B. von Issendorff, H. Haberland, P. Labastie, Phys. Rev. Lett. 87 (2001) 203402. [18] P. Parneix, Ph. Br
echignac, J. Chem. Phys. 118 (2003) 8234. [19] M.F. Jarrold, in: H. Haberland (Ed.), Clusters of Atoms and Molecules I, Springer, Berlin, 1991. [20] F. Wang, D.P. Landau, Phys. Rev. Lett. 86 (2001) 2050.

[21] F. Calvo, Mol. Phys. 100 (2002) 3421. [22] R.H. Swendsen, J.-S. Wang, Phys. Rev. Lett. 57 (1986) 2607; G.J. Geyer (Ed.), Computing Science and Statistics: Proceedings of the 23rd Symposium on the Interface, American Statistical Association, New York, 1991. [23] A.M. Ferrenberg, R.H. Swendsen, Phys. Rev. Lett. 61 (1988) 2635. [24] D.J. Wales, Mol. Phys. 78 (1993) 151. [25] V. Weisskopf, Phys. Rev. 52 (1937) 295. [26] P. Parneix, Eur. Phys. J. D 23 (2003) 375.