TEMPERATURE INDUCED PHASE TRANSFORMATIONS OF MOLECULAR NANOCLUSTERS

A. Proykova^{1},
S. Pisov^{1}, R. Radev^{1}, P. Mihailov^{1},
I. Daykov^{2}, and R.S.Berry^{3}

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

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

^{3}The 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 TeF_{6}
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 TeF_{6
}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 TeF_{6}
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 TeF_{6 }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 (TeF_{6})
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 TeF_{6} -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
TeF_{6} 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.1*e* 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]:

(2)

where ** **(

**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 *k _{ii}* 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

**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].

**Conclusion**

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.

**References**

Berry RS, Beck, Davis HL, Jellinek J. in

*Evolution of Size Effects in Chemical Dynamics*, Part 2, (New-York, 1988); p.75Bartell LS. J. Phys. Chem. (1992);

**96**:108Proykova A, Berry RS. Z. Phys D. (1997);

**40**:215Proykova A, Radev R, Li FY, Berry RS. J. Chem. Phys. (1999);

**110**:3887Bartell LS, Valente EJ, Caillat J. J. Phys. Chem. (1987);

**91**:2498.Allen MP, Tildesley DJ.:

*Computer Simulation of Liquids*( Clarendon Press, Oxford 1994).Stillinger FH, Weber TA.

*Kinam***3**, 159 (1981).Proykova A, Pisov S, Berry RS. J. Chem.Phys. (2001);

**115**:8583Binder K.

*Monte Carlo Methods in statistical physics*. (Springer - Berlin 1979).Binney JJ, Dowrick NJ, Fisher AJ, Newman MEJ.

*The theory of Critical Phenomena*, Oxford Science Publications, (Clarendon Press - Oxford 1992).Pisov S, Proykova A.

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

*Meeting in Physics at University of Sofia,*ed. A. Proykova,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*.Mihailov P.

*MS Diploma thesis*(1998).Frantz DD, Freeman DL, Doll JD. J. Chem. Phys. (1990);

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

**97**:5713Radev R.

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

*Proceedings of the Conference on Computational Physics*(CCP01*), Aachen,*5-9 Sep., Germany (2001).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 TeF_{6
}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 TeF_{6} cluster at T= 157 K
(liquid) and T= 76 K (two solid structures can be distinguished).

**Fig. 9 **The**
**caloric curve of a 59 TeF_{6} 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.1 |

Fig.2 |

Fig.3 |

Fig.4 |

Fig.5 Fig.6

Fig.7

Fig.8

Fig.9