A. Proykova1, S. Pisov1, R. Radev1, P. Mihailov1, I. Daykov2, and R.S.Berry3

1University of Sofia, Faculty of Physics, 5 J. Bourchier Blvd., Sofia-1126, Bulgaria

2Cornell University, Department of Physics, 117 Clark Hall, Ithaca, NY 14850, USA

3The University of Chicago, Department of Chemistry, 5735 S. Ellis Ave., Chicago IL, 60637, USA

Abstract: We report results from canonical Molecular Dynamics (MD) and Jump-Walking Monte Carlo (MC) simulations, which show two sequential solid-solid transformations of nanosized (~ 2nm) free TeF6 clusters on cooling below their solidification point. MC and MD simulations define the same temperature interval (76±4) K for the discontinuous (first order) transformation of 59-molecule TeF6 clusters and (29±5) K for the continuous transformation. We conclude that the ergodic requirements can be satisfied in solid-solid transformations of these specific clusters. The discontinuous transition is characterized with dynamically coexisting phases in a single cluster observed in MD simulations. The coexisting phases are stabilized by kinetics of small systems, although they might be thermodynamically unstable in the bulk of the same material.

Key Words: nano-sized clusters, phase changes, molecular dynamics, jump-walking Monte Carlo simulations

1. Introduction

Clusters have become increasingly popular systems to investigate the kind and mechanism of phase changes, [1-4]. Transitions between different phases are not only commonly encountered in nature but they play a central role in many technological processes in industry, and are also of fundamental interest for theoretical physics in understanding the behavior of large number of interacting particles.

Previously, we have used microcanonical Molecular Dynamics (MD) method to study temperature driven phase-like transformations of TeF6 molecule clusters of various sizes [3,4]. Liquid-like and solid-like phases can be distinguished by the value of the Lindemann index, [1]. In bulk, < 0.1 indicates a solid state. In a free cluster the surface plays a destabilizing role and a cluster is solid if< 0.08, otherwise it is melted. The different solid structures are distinguished on the basis of their radial distributions of the molecular mass g(r), [6].

It has been shown that varieties of structures appear below the freezing point, which depends on the cluster size: the smaller the cluster the lower the temperature [2-5]. This abundance of structures reflects the complex intermolecular interaction that depends on the space orientation of the molecules. The potential energy surface (PES) contains many minima with the same depth, which corresponds to many isomers found by quenching of the cluster along its trajectory. Quenching is a means to associate arbitrary points in the phase space with the local minima in whose wells they lie on the potential energy surface, [7].

In the present study we perform both MD and Monte Carlo (MC) canonical simulations of temperature driven phase changes of free molecular clusters containing 59 TeF6 molecules. This number is the smallest one, which ensures one volume cell having all the neighbors required by the symmetry of the body centered cubic lattice observed in the experiments of cluster solidification [5]. The aim of the present study is twofold. Firstly, we check if the prediction of the two methods coincide within the calculation errors, e.g. to check the validity of the ergodic hypothesis in small systems that have rugged potential energy surfaces. Secondly, we investigate whether distinct coexisting structures can be dynamically observed in one and the same small cluster while in bulk they are observable in ensembles. The second subject is very important for understanding the properties of cluster-assembled materials (nanopowders, polymers) used in the modern nanotechnology.

Clusters made of tellurium hexafluoride (TeF6) molecules have rugged PES [8]. Such a PES contains many minima separated by high barriers that hamper the expected evolution of the system at a given temperature. The freezing kinetics of TeF6 -molecule clusters with different sizes has been examined in some detail, [2-5]. The main difficulty in the analysis of the system properties is the so-called trapping. The cluster trajectory in the phase space could be confined in a local minimum for a long time, thousands of characteristic vibration times, , defined below. If this happens, one gets biased values of the time-averaged quantities such as internal energy or entropy. Increasing the time of computations is a possible solution to this problem, however, one should not forget that very long runs violate the energy conservation even in the most efficient integration schemes, like simplex method [8]. 'Violation' means that if the same computation is performed with t replaced with (-t) the trajectory will not return to its initial point in the phase space.

State averaging of the internal energy, for example, obtained by Monte Carlo methods can give a complementary information about the system [9]. If the system is ergodic, or at least approximately ergodic [10], the time- and state-averages should be the same in the limits of the computational errors. The system is considered to be ergodic if it passes, with equal probability, through all states available to the system with a given energy. Ergodicity requires that there are no other constants of motion besides the energy. Molecular clusters exhibiting order-disorder phase changes do not conserve the angular momentum of the molecules due to alignment of molecular axes in the ordered phase [3]. Thus we can assume that the system is naturally ergodic. However, the computer simulations are limited in time and we do not know in advance if the runs performed have been long enough to retrieve the system properties correctly. The predictions of MD and MC might be quite different. That is why it is necessary to check the validity of the ergodic hypothesis in the simulations of temperature driven phase transitions of molecular nanoclusters by implementing both dynamic and static methods for the same system governed by a specific potential.

In section 2 we consider a molecular substance whose ingredients interact with an orientation dependent potential that generates a rugged potential surface. Section 3 presents results from canonical MD computations. Section 4 briefly describes the Jump-Walking Monte Carlo technique used in our simulations. Concluding remarks complete the article.

2. Molecular Clusters with Rugged Potential Energy Surfaces

The TeF6 molecules can be considered as rigid octahedrons: the positively charged tellurium atom is in the center and the six, negatively charged fluorine atoms are symmetrically located at a constant distance, the bond length of 1.814 Å. The fluorine charge of 0.1e has been computed using the Linear Combination of Atomic Orbitals (LCAO) approximation with planar basic functions [3]. A Gaussian basis gives a charge of 0.25 e [11]. The increased charge changes the PES topography: the wells become deeper and the transitions between the solid structures are suppressed as our study indicates [11].

To account for its orientation dependence the total potential is presented as a sum over atom-atom potentials of a Lennard-Jones and a Coulomb type, [3]:

U(r)=, (1)

where is the distance between the i-th and j-th atom in different molecules. There is no intramolecular interaction in the rigid molecules. The indices , denote either fluorine or a tellurium atom. The parameters are = 4.29 Å, = 3.55 Å, = 2.94 Å, F-F = 0.2049 kJ/mol, F-Te = 1.376 kJ/mol, Te-Te = 0.531 kJ/mol. The parameter values have been fitted to match the diffraction patterns obtained in the experiments, [2]. Due to long-range interactions and the small cluster size no cutoffs of the potential have been applied in the calculations.

The plot in Fig.1 shows the potential energy surface of a cluster containing 3 molecules arranged as follows. Two of the molecules are fixed at a distance of 6.18 Å, corresponding to the potential minimum of two aligned molecules. The third molecule, aligned to the others, probes the potential surface at different distances from the cluster center of masses. The picture is inverted, so the potential minima are seen as peaks. The highest peak corresponds to the most stable structure. However, there are several peaks with the same heights, which means that there are several isomers of the studied cluster. Hence, a cluster formed as a specific isomer could not transform to another at a specific temperature in a short time interval. If the cluster is heated, just a little bit, and after that it is slowly cooled, then a transformation to the most stable structure is possible.

The temperature behavior of the clusters has been determined by solving the Hamilton equations of motion written in the Nose’s form [12]:


where (r,p) are the coordinates and momenta of the particles, i = 1,…., N, N is the number of the molecules in the cluster, g is the number of degrees of freedom, kB is the Boltzman constant, T is the temperature, M has the meaning of a bath size, which governs the heat transfer between the system and the bath. The 'friction' controls the temperature of the system . The bath size M depends both on the system size and the interacting potential U. For the potential given with Eq.1 and 59 TeF6 molecules, a proper bath size is M = 200. The right choice of the parameter M is very important for the temperature of the cluster. If it is too small the new degree of freedom decouples from the cluster and has its own dynamics, resulting in a double-peak distribution of the kinetic energy (like a harmonic oscillator) rather than the desired Gaussian distribution. If it is too big the heat bath is very heavy and one has to wait too long for establishing thermal equilibrium. A detailed study of the bathsize for various sizes of TeF6 and SF6 molecule clusters is published elsewhere [12].

3. Results from Molecular Dynamics Computations

Variety of structures originates from the orientation dependence of the intermolecular potential, Eq.1, when the temperature of the cluster is decreased below the freezing point: first, a body-centered cubic structure is observed in real and computer experiments, [2,4,5].

The molecular axes of symmetry are orientationally disordered as it is seen in Fig.2. Further cooling of the clusters causes solid-solid transformation to a monoclinic or an orthorhombic lattice depending on the cooling rate. The molecules are completely orientationally ordered below 20K, as it is well seen in the quenched structure, Fig.3.

The quenched structure at the transition temperature, which depends on the actual cluster size, reveals two coexisting structures, Fig.4. We have observed coexisting structures in many clusters at a properly chosen cooling (heating) rate only. The change of the temperature in the vicinity of phase-like transformations should be about 1K for the 59 molecule clusters and much smaller for bigger clusters: 0.1K for 137 molecule clusters.

In fact, the coexistence is not observed at a unique temperature but in a temperature interval, which becomes narrower as the cluster grows larger, [12]. Clusters containing 59 molecules transform from a body centered, orientationally disordered phase to a partially ordered monoclinic phase at (76 ±5) K. This transformation is detected at (89±4) for clusters of 89 molecules and at (91.1±0.9) K for 137-molecule clusters. The temperature interval is determined by successive heating of clusters that have been cooled: a hysteresis is observed in all simulated clusters (about 100) of the same size [4]. By averaging over hysteresis widths we obtain the temperature interval. Appearance of a finite temperature interval of the phase transformation is a phenomenon that distinguishes the finite systems from their bulk counterpart that transforms at a unique temperature (233K for the bulk tellurium hexafluoride). The phase coexistence of two distinct structures indicates that the phase transformation is of first order (discontinuous). The free energies of the two phases are similar enough for them to coexist in observable amounts in the transition temperature interval while their entropies and internal energies are different. The plot of the internal energy (see Figs.3 and 4 in ref. [13]) reveals two distinct states visited by the system during its evolution at the transition temperature.

Coexisting phases are dynamically observed in simulation runs much longer than the characteristic vibrational time , where kii are the diagonal elements of the Hessian matrix, . The rotation (libration) characteristic time is shorter than the intermolecular vibration period, which is about 2.5 ps for the clusters of interest. As it has been mentioned, we need at least 0.5 ns in order to observe coexisting phases in 59 molecule clusters. For larger clusters, the observation time should be correspondingly increased [12]. As example, the time needed is 5 ns for a 137 molecule cluster.

4. Monte Carlo Techniques for Sampling Rugged Potential Energy Surfaces

The Monte Carlo simulations give complementary information about the studied system by generating a Markov sequence of configurations. These are used to obtain the statistical (ensemble) average of the measured quantity [9]. In a molecular dynamic simulation we compute a long trajectory in the phase space over which we calculate the time averages of the same quantities. Under the conditions of the ergodic hypothesis these two should be the same, e.g. <A>{time} = <A>{configurations}.

At low temperatures, the conventional Metropolis Monte Carlo methods [14] might not properly sample configurations of systems described with a rugged PES. We have found that extremely long runs (millions of Monte Carlo steps) sometimes do not let the system to equilibrate to the lowest-lying state, which we know from the molecular dynamics simulations.

An improvement is achieved by the jump-walking technique for spanning the configuration space, [15, 16]. This technique mixes Boltzman weighted states, corresponding to low and high temperatures. We have developed [17] a parallel J-Walking code that enables us to carry out a Monte Carlo simulation efficiently in a multi-processor-computing environment and have applied it in the study of the thermodynamic properties of 59-molecule clusters.

A schematic representation of the Jump-Walking (J-W) procedure is given in Fig.5. Each square represents a Metropolis Monte Carlo simulation at a particular temperature. The boxes in the right-hand side denote the array of previous configurations stored in the memory to avoid correlation between the lower and the higher temperature. Such a procedure is equivalent to the following physical situation: a cluster at a given temperature could contain grains, which do not let it transform into the most stable structure. In order to get rid of them, we shortly heat the cluster to let the domains melt their boundaries. The algorithm does the same by mixing different states into the same simulation. The essence of mixing is that we increase the probability of population of states, which are hardly accessible for the cluster before heating. Such kind of techniques are needed to explore rugged potential surfaces that usually trap the cluster in a specific well, normally corresponding to a glassy structure. Recently, we have found that the parallel tempering technique is even better than the J-Walking Monte Carlo for achieving equilibrium of the system, [18].

In order to obtain reliable results we need a great many configurations, at least several millions. A speed-up of computations is achieved by using parallel computers. In Fig.6 we give an example of the parallel decomposition of simultaneous simulations of 4 clusters at four temperatures distributed between 8 processors. Each processor computes one of the four multistage J-walking and exchanges configurations and energy with the others.

The parallel MC code has been implemented to simulate the temperature behavior of four 59-molecule clusters at 32 different temperatures distributed between 64 processors of the parallel machines SUN3500 and CRAY-T3E.

We use the Lindemann index to distinguish between liquid-like and solid-like forms and as well, between two different solid structures. In Fig.7 we plot the Lindemann index as a function of the temperature. Curves corresponding to 4 different clusters are given. Melting of these clusters starts at ~100 K. Below this temperature, (T) changes its slope at ~80K (the first solid-solid transformation) indicating a transition between two solid structures because its values are much lower than 0.08. Another indicator of the cluster states (liquid, solid) is the radial distribution of the molecular centers of masses, g(r). The liquids are characterized with g(r) having regularly located peaks [6] as it is seen in Fig.8 for the case of T= 157 K, much above the melting point. The difference between a liquid-like cluster and bulk liquid is the decreasing peak heights due to the cluster finite size. At T = 76 K, the clusters are solid and the curve indicates a mixture of bcc and monoclinic structures. Below this temperature, at ~ 30 K we detect another transformation in all simulations (Molecular Dynamics, conventional Metropolis and J-Walking Monte Carlo). In Fig.9 we present the caloric curve of clusters simulated with the conventional Metropolis MC method. The temperature is well below the discontinuous solid-solid transition at ~80 K. In the region 25 -35 K the curve changes its slope and the fluctuations of the mean potential energy increase. Another transition occurs - from partially aligned to completely aligned molecular axes, as it is seen in Fig. 3. In bulk we would observe a smooth change of the caloric curve in (25-35) K region, as it is typical for the continuous transformations [10]. The temperature region with increased fluctuations of the internal energy is taken as the temperature interval of the transformation. It is quite difficult to determine the kind of transitions in finite systems only from simulations. To confirm that the second transition is continuous we have applied a point-group approach to the solid-solid transformations of molecular clusters [19].


Both Molecular Dynamics and Monte Carlo methods show two solid-solid transformations of 59 TeF6-molecule clusters below their freezing point: discontinuous transition in (76±4)K and continuous transition in (29±5)K. The same transition temperature interval for the cluster size studied (~2.4 nm) is obtained if long MD runs are carried out (~ 1 ns) and several millions of configurations are used in MC. The interval width is determined by averaging over many simulations of different clusters with the same size. We conclude that the ergodic requirements can be satisfied in solid-solid transformations of nanoclusters under those special limitations.

J-Walking Monte Carlo techniques allows the clusters to reach equilibrium much faster than the conventional Metropolis Monte Carlo method and decreases the computational efforts.

Acknowledgements: The work has been completed under the NATO Collaborative Linkage Grant Award Number: SA(PST.CLG.976363)5437; the TRACS EU program (Edinburgh) supported R.Radev to develop the parallel Monte Carlo J-Walking code; support from the Cornell University (I. Daykov) and from the University of Sofia (S. Pisov, University of Sofia Scientific Fund No. 3270/2001) is greatly acknowledged.


  1. Berry RS, Beck, Davis HL, Jellinek J. in Evolution of Size Effects in Chemical Dynamics, Part 2, (New-York, 1988); p.75

  2. Bartell LS. J. Phys. Chem. (1992); 96:108

  3. Proykova A, Berry RS. Z. Phys D. (1997); 40:215

  4. Proykova A, Radev R, Li FY, Berry RS. J. Chem. Phys. (1999); 110:3887

  5. Bartell LS, Valente EJ, Caillat J. J. Phys. Chem. (1987); 91:2498.

  6. Allen MP, Tildesley DJ.: Computer Simulation of Liquids( Clarendon Press, Oxford 1994).

  7. Stillinger FH, Weber TA. Kinam 3, 159 (1981).

  8. Proykova A, Pisov S, Berry RS. J. Chem.Phys. (2001); 115:8583

  9. Binder K. Monte Carlo Methods in statistical physics. (Springer - Berlin 1979).

  10. Binney JJ, Dowrick NJ, Fisher AJ, Newman MEJ. The theory of Critical Phenomena, Oxford Science Publications, (Clarendon Press - Oxford 1992).

  11. Pisov S, Proykova A. Proceedings of the Conference on Computational Physics (CCP01), Aachen, Germany (2001); B21. (to appear in Comp. Phys. Com.)

  12. Daykov I, Proykova A. Meeting in Physics at University of Sofia, ed. A. Proykova, v.2 p.31(Heron Press Science Series - Sofia 2001).

  13. Proykova A. Proceedings of the Second workshop on Nanoscience and Nanotechnology, eds. E. Balabanova and I. Dragieva, ISBN 954-580-097-6 (Heron Press – Sofia 2001); Proykova A. submitted to Vacuum.

  14. Mihailov P. MS Diploma thesis (1998).

  15. Frantz DD, Freeman DL, Doll JD. J. Chem. Phys. (1990); 93:2769

  16. Frantz DD, Freeman DL, Doll JD. J. Chem. Phys. (1992); 97:5713

  17. Radev R. Parallel Monte Carlo simulation of critical behavior of finite systems, Report from the TRACS program, Edinburgh-UK (May-June 2000).

  18. Radev R, Proykova A. Proceedings of the Conference on Computational Physics (CCP01), Aachen, 5-9 Sep., Germany (2001). (to appear in Com.Phys.Com.)

  19. Proykova A, Nikolova D, Berry RS. (to appear in Phys. Rev. B); preprint http://xxx.lanl.gov/abs/physics/9911049

Figure Captions

Fig.1 Potential energy surface of a cluster containing 3 TeF6 molecules.

Fig.2 Orientationally disordered body centered cubic lattice seen below the freezing point.

Fig.3 A stable orientationally ordered structure is detected below the second solid-solid transformation in (29±5) K.

Fig.4 Coexisting structures observed in (76±4)K signal discontinuous transition.

Fig.5 The diagram shows the J-W process between two temperatures.

Fig.6 Parallel decomposition of data and calculation among the processors in J-walking algorithm.

Fig.7 The Lindemann index shows that the cluster is solid below 100K and it transforms into another solid structure at ~80 K.

Fig. 8 Radial distribution of a 59 TeF6 cluster at T= 157 K (liquid) and T= 76 K (two solid structures can be distinguished).

Fig. 9 The caloric curve of a 59 TeF6 cluster shows a solid-solid transformation in the vicinity of 30 K. The result has been obtained with the help of the conventional Metropolis Monte Carlo method.





Fig.5 Fig.6