Phase changes in 38 atom Lennard-Jones clusters. I: A parallel tempering study in the canonical ensemble

The heat capacity and isomer distributions of the 38 atom Lennard-Jones cluster have been calculated in the canonical ensemble using parallel tempering Monte Carlo methods. A distinct region of temperature is identified that corresponds to equilibrium between the global minimum structure and the icosahedral basin of structures. This region of temperatures occurs below the melting peak of the heat capacity and is accompanied by a peak in the derivative of the heat capacity with temperature. Parallel tempering is shown to introduce correlations between results at different temperatures. A discussion is given that compares parallel tempering with other related approaches that ensure ergodic simulations.

which, unlike the case for most small Lennard-Jones clusters, is not based on an icosahedral core, but rather is a symmetric truncated octahedron. The vertices defined by the surface atoms of LJ 38 have a morphology identical to the first Brillouin zone of a face centered cubic lattice, 10 and the high symmetry of the cluster may account for its stability. It is interesting to note that recent experimental studies 11 of nickel clusters using nitrogen uptake measurements have found the global minimum of Ni 38 to be a truncated octahedron as well.
The basin of energy minima about the global minimum of LJ 38 is narrow compared to the basin about the next highest energy isomer which does have an icosahedral core. The difference in energy between the global minimum and the lowest minimum in the icosahedral basin is only 0.38% of the energy of the global minimum. 5 Characteristic of some thermodynamic properties of small clusters are ranges of temperature over which these properties change rapidly in a fashion reminiscent of the divergent behavior known to occur in bulk phase transitions at a single temperature. The rapid changes in such thermodynamic properties for clusters are not divergent and occur over a range of temperatures owing to the finite sizes of the systems. In accord with the usage introduced by Berry, Beck, Davis and Jellinek 12 we refer to the temperature ranges where rapid changes occur as "phase change" regions, rather than using the term "phase transition," that is reserved for systems at the thermodynamic limit. As an example LJ 55 displays a heat capacity anomaly over a range in temperatures often associated with what has been termed "cluster melting." 13 Molecular dynamics and microcanonical simulations performed at kinetic temperatures in the melting region of LJ 55 exhibit van der Waals type loops in the caloric curves and coexistence between solid-like and liquid-like forms.
In recent studies, Doye, Wales and Miller 14 and Miller, Doye and Wales 15 have examined the phase change behavior of LJ 38 . These authors have calculated the heat capacity and isomer distributions as a function of temperature using the superposition method. 16,17 In the superposition method the microcanonical density of states is calculated for each potential minimum, and the total density of states is then constructed by summation with respect to each local density of states. Because it is not possible to find all potential minima for a system as complex as LJ 38 , the summation is augmented with factors that represent the effective weights of the potential minima that are included in the sum. The superposition method has also been improved to account for anharmonicities and stationary points. 17 For LJ 38 Doye et al. 14 have identified two phase change regions. The first, accompanied by a heat capacity maximum, is associated with a solid-to-solid phase change between the truncated octahedral basin and the icosahedral basin. A higher temperature heat capacity anomaly represents the solid-liquid coexistence region, similar to that found in other cluster systems.
The heat capacity anomaly associated with the melting transition in LJ 38 is steeper and more pronounced than the heat capacity peak that Doye et al. 14 have associated with the solidsolid transition. Because the weights that enter in the sum to construct the microcanonical density of states are estimated, it is important to confirm the findings of Doye et al. 14 by detailed numerical simulation. Such simulations are a goal of the current work and its companion paper. As is found in Section III, the simulations provide a heat capacity curve for LJ 38 that has some qualitative differences with the curve reported by Doye et al. 14 Owing to the complex structure of the potential surface of LJ 38 , the system represents a particularly challenging case for simulation. It is well known that simulations of systems having more than one important region of space separated by significant energy barriers can be difficult. The difficulties are particularly severe if any of the regions are either narrow or reachable only via narrow channels. The narrow basin about the global minimum makes simulations of LJ 38 especially difficult. There are several methods that have been developed that can prove to be useful in overcoming such ergodicity difficulties in simulations. Many of these methods use information about the underlying potential surface generated from simulations on the system using parameters where the various regions of configuration space are well-connected. One of the earliest of these methods is J-walking 18 where information about the potential surface is obtained from simulations at high temperatures, and the information is passed to low temperature walks by jumping periodically to the high temperature walk.
Closely allied with J-walking is the parallel tempering method [19][20][21][22][23] where configurations are exchanged between walkers running at two differing temperatures. A related approach, 24 similar in spirit to J-walking, uses Tsallis distributions that are sufficiently broad to cover much of configuration space. Another recent addition 25 to these methods is the use of multicanonical distributions 26 in the jumping process. Multicanonical walks are performed using the entropy of the system, and multicanonical distributions are nearly independent of the energy thereby allowing easy transitions between energy basins. As we discuss in the current work, we have found the parallel tempering method to be most useful in the context of simulations of LJ 38 . A comparative discussion of some of the methods outlined above is given later in this paper.
In the current work we apply parallel tempering to the calculation of the thermodynamic properties of LJ 38 in the canonical ensemble. In the paper that follows 27 we again use parallel tempering to study LJ 38 , but using molecular dynamics methods along with microcanonical Monte Carlo simulations. Our goals are to understand better this complex system and to determine the best simulation method for systems of comparable complexity. The contents of the remainder of this first paper are as follows. In Section II we discuss the methods used with particular emphasis on the parallel tempering approach and its relation to the J-walking method. In Section III we present the results including the heat capacity as a function of temperature and identify the phase change behaviors of LJ 38 . In Section IV we present our conclusions and describe our experiences with alternatives to parallel tempering for insuring ergodicity.

II. METHOD
For canonical simulations we model a cluster with N atoms by the standard Lennard-Jones potential augmented by a constraining potential U c used to define the cluster where σ and ε are respectively the standard Lennard-Jones length and energy parameters, and r ij is the distance between particles i and j. The constraining potential is necessary because clusters at defined temperatures have finite vapor pressures, and the evaporation events can make the association of any atom with the cluster ambiguous. For classical Monte Carlo simulations, a perfectly reflecting constraining potential is most convenient with where r cm is the center of mass of the cluster, and we call R c the constraining radius.
Thermodynamic properties of the system are calculated with Monte Carlo methods using the parallel tempering technique. [19][20][21][22][23] To understand the application of the parallel tempering method and to understand the comparison of parallel tempering with other related methods, it is useful to review the basic principles of Monte Carlo simulations.
In the canonical ensemble the goal is the calculation of canonical expectation values. For example, the average potential energy is expressed where β = 1/k B T with T the temperature and k B the Boltzmann constant. In Monte Carlo simulations such canonical averages are determined by executing a random walk in configuration space so that the walker visits points in space with a probability proportional to the canonical density ρ(r) = Z −1 exp[−βU(r)], where Z is the configurational integral that normalizes the density. After generating M such configurations in a random walk, the expectation value of the potential energy is approximated by The approximate expectation value U M becomes exact in the limit that M → ∞.
A sufficiency condition for the random walk to visit configuration space with a probability proportional to the density ρ(r) is the detailed balance condition 28,29 where r o and r n represent two configurations of the system and K(r o → r n ) is the conditional probability that if the system is at configuration r o it makes a transition to r n . In many Monte Carlo approaches, the conditional probability is not known and is replaced by the expression where T (r o → r n ) is called the trial probability and acc(r o → r n ) is an acceptance probability constructed to ensure K(r o → r n ) satisfies the detailed balance condition. The trial probability can be any normalized density function chosen for convenience. A common choice for the acceptance probability is given by 28,29 acc(r o → r n ) = min 1, Metropolis method rigorously guarantees a random walk visits configuration space proportional to a given density function asymptotically in the limit of an infinite number of steps.
In practice when configuration space is divided into important regions separated by significant energy barriers, a low temperature finite Metropolis walk can have prohibitively long equilibration times.
Such problems in attaining ergodicity in the walk do not occur at temperatures sufficiently high that the system has significant probability of finding itself in the barrier regions.
In both the J-walking and parallel tempering methods, information obtained from an ergodic Metropolis walk at high temperatures is passed to a low temperature walker periodically to enable the low temperature walker to overcome the barriers between separated regions. In the J-walking method 18 the trial probability at inverse temperature β is taken to be a high temperature Boltzmann distribution where β J represents the jumping temperature that is sufficiently high that a Metropolis walk can be assumed to be ergodic. Introduction of Eq.(9) into Eq.(8) results in the acceptance In practice at inverse temperature β the trial moves are taken from the Metropolis distribution about 90% of the time with jumps attempted using Eq.(9) about 10% of the time.
The jumping configurations are generated with a Metropolis walk at inverse temperature the probability of ever choosing the same configuration more than once is small. In this method detailed balance is strictly satisfied only in the limit that the external distribution is of infinite size. In the second method, often called parallel J-walking, 31,32 the walks at each temperature are made in tandem on a parallel machine. Many processors, randomly initialized, are assigned to the jumping temperature, and each processor at the jumping temperature is used to donate a high temperature configuration to the low temperature walk sufficiently infrequently that the correlations in the Metropolis walk at inverse temperature β J are broken. In this parallel method, configurations are never reused, but the acceptance criterion [Eq. (10)] is strictly valid only in the limit of an infinite set of processors at inverse temperature β J . In practice both serial and parallel J-walking work well for many applications with finite external distributions or with a finite set of processors. 18,[31][32][33][34][35][36][37][38] In parallel tempering [19][20][21][22][23] configurations from a high temperature walk are also used to make a low temperature walk ergodic. In contrast to J-walking rather than the high temperature walk feeding configurations to the low temperature walk, the high and low temperature walkers exchange configurations. By exchanging configurations detailed balance is satisfied, once the Metropolis walks at the two temperatures are sufficiently long to be in the asymptotic region. To verify detailed balance is satisfied by the parallel tempering procedure we let be the joint density that the low temperature walker is at configuration r and the high temperature walker is at configuration r ′ . When configurations between the two walkers are exchanged, the detailed balance condition is By solving for the ratio of the conditional transition probabilities it is evident that if exchanges are accepted with the same probability as the acceptance criterion used in J-walking [see Eq.(10)], detailed balance is satisfied.
Although the basic notions used by both J-walking and parallel tempering are similar, the organization of a parallel tempering calculation can be significantly simpler than the organization of a J-walking calculation. In parallel tempering no external distributions are required nor are multiple processors required at any temperature. Parallel tempering can be organized in the same simple way that serial tandem J-walking is organized as discussed in the original J-walking reference. 18 Unlike serial tandem J-walking where detailed balance can be attained only asymptotically, parallel tempering satisfies detailed balance directly.
For a problem as difficult as LJ 38 where very long simulations are required, the huge external distributions needed in serial J-walking, or the large set of jumping processors needed in parallel J-walking, make the method prohibitive. As discussed in Section III, parallel tempering can be executed for arbitrarily long simulations making the method suitable at least for LJ 38 .
In the current calculation parallel tempering is used not just to simulate the system at some low temperature using high temperature information, but simulations are performed for a series of temperatures. As is the case for J-walking 18 and as discussed elsewhere for parallel tempering, 22 the gaps between adjacent temperatures cannot be chosen arbitrarily.
Temperature gaps must be chosen so that exchanges are accepted with sufficient frequency.
If the temperature gap is too large, the configurations important at the two exchanging temperatures can be sufficiently dissimilar that no exchanges are ever accepted. Preliminary calculations must be performed to explore the temperature differences needed for acceptable exchange probabilities. In practice we have found at least 10% of attempted exchanges need to be accepted for the parallel tempering procedure to be useful. In general the temperature gaps must be decreased near phase change regions or when the temperature becomes low.
By exchanging configurations between temperatures, correlations are introduced at different temperature points. For example, the average heat capacities at two temperatures may rise or fall together as each value fluctuates statistically. In some cases the values of the heat capacities or other properties at two temperatures can be anti-correlated. The magnitude of these correlations between temperatures are measured and discussed in Section III. As discussed in Section III the correlations between differing temperatures imply that the statistical fluctuations must be sufficiently low to ensure any features observed in a calculation as a function of temperature are meaningful.

III. RESULTS
Forty distinct temperatures have been used in the parallel tempering simulations of LJ 38 In an attempt to minimize the correlations in the data at differing temperatures, an exchange strategy has been used that includes exchanges between several temperatures. To understand this strategy, we let the set of temperatures be put into an array. One-half of the exchanges have been attempted between adjacent temperatures in the array, one-fourth have been attempted between next near neighboring temperatures, one-eighth between every third temperature, one-sixteenth between every fourth temperature and one-thirty second between every fifth temperature in the array. We have truncated this procedure at fifth near neighboring temperatures, because exchanges between temperatures differing by more than fifth neighbors are accepted with frequencies of less than ten per cent. The data presented in this work have been generated using the procedure outlined above. In retrospect, we have is displayed in the upper panel of The small low temperature maximum in (∂C V /∂T ) V occurs within the slope change region.
To interpret the configurations associated with the various regions of the heat capacity, we use an order parameter nearly identical to the order parameter introduced by Steinhardt, Nelson and Ronchetti 39 to distinguish face centered cubic from icosahedral structures in liquids and glasses. The order parameter has been used by Doye et al. 5 to monitor phase changes in LJ 38 . The order parameter Q 4 is defined by the equation where To understand the parameters used in Eq. (17), it is helpful to explain how Q 4,m is evaluated.
The center of mass of the full 38 atom cluster is located and the atom closest to the center of mass is then identified. The atom closest to the center of mass plus the 12 nearest neighbors of that atom define a "core" cluster of the 38 atom cluster. The center of mass of the core cluster is then calculated. The summation in Eq. (17) is performed over all vectors that point from the center of mass of the core cluster to all N b bonds formed from the 13 atoms of the core cluster. A bond is assumed to be formed between two atoms of the core cluster if their internuclear separation r ij is less than a cut-off parameter r b , taken to be r b = 1.39σ in this work. In Eq.(17) θ ij and φ ij are respectively the polar and azimuthal angles of the vector that points from the center of mass of the core cluster to the center of each bond, and Y 4,m (θ, φ) is a spherical harmonic. To verify that the optimal value of Q 4 is obtained, the procedure is repeated by choosing the second closest atom to the center of mass of the whole cluster to define the core cluster. The value of Q 4 obtained from this second core cluster is compared with that obtained from the first core cluster, and the smallest resulting value of Q 4 is taken to be the value of Q 4 for the entire cluster.
In the work of Steinhardt et al. 39   as a function of temperature and order parameter. A projection of this surface onto two dimensions is given in Fig. 5. The probability density in Fig. 5 is represented by the shading so that the brighter the area the greater the probability. The horizontal white lines in Fig. represents the temperature at which the slope of the heat capacity first changes rapidly, the middle temperature horizontal line represents the lowest temperature of the melting peak and the highest temperature horizontal line represents the end of the melting region. An additional representation of the data is given in Fig. 6, where the probability of observing particular values of Q 4 is given as a function of Q 4 at a fixed temperature of 0.14ε/k B . In  There is significant density for both kinds of structures in the region between the lowest two parallel lines that define the region with the slope change. Additionally, both icosahedral structures and truncated octahedral structures begin to be in equilibrium with each other at the beginning of the slope change region. This equilibrium continues to temperatures above the melting region.
Another identification of the slope change region with a transition between truncated octahedral and icosahedral forms can be made by defining P R (T, R)dR to be the probability that an atom in the cluster is found at location R to R + dR from the center of mass of the cluster at temperature T . A projection of P R (T, R) onto the R and T plane is depicted in does not match any of the vertical lines shown, but corresponds to atoms in the third lowest energy isomer, which like the second lowest energy isomer, comes from the icosahedral basin.
The equilibrium between the icosahedral and truncated octahedral structures observed in Fig. 7 matches the regions of temperature observed in Fig. 5.
We have mentioned previously that parallel tempering introduces correlations in the data accumulated at different temperatures, and it is important to ensure the statistical errors are sufficiently small that observed features are real and not artifacts of the correlations.
To measure these correlations we define a cross temperature correlation function for some temperature dependent property g by A projection of γ(T 1 , T 2 ) when g = C V is given in Fig. 8. In Fig. 8  heat capacity maximum associated with this transition, but rather a region with a change in the slope of the heat capacity as a function of temperature.
We have found parallel tempering to be successful with this system, and have noted correlations in our data at different temperatures when the parallel tempering method is used. These correlations imply the need to perform long simulations so that the statistical errors are sufficiently small that the correlations do not introduce artificial conclusions.
We believe that the methods used in this work could be applied to a variety of other