A j-Walking Algorithm for Microcanonical Simulations: Applications to Lennard-Jones Clusters

The j-walking method, previously developed to solve quasiergodicity problems in canonical simulations, is extended to simulations in the microcanonical ensemble. The implementation of the method in the microcanonical ensemble parallels that in the canonical case. Applications are presented in the microcanonical ensemble to cluster melting phenomena for Lennard-Jones clusters containing 7 and 13 particles. Signiﬁcant difﬁculties are encountered in achieving ergodicity when Metropolis Monte Carlo methods are applied to these systems, and the difﬁculties are removed by the j-walking method. © 1998 American Institute of Physics. (cid:64) S0021-9606 (cid:126) 98 (cid:33) 50129-8 (cid:35)


I. INTRODUCTION
The j-walking method 1,2 has proved to be useful in overcoming quasiergodicity problems for Metropolis Monte Carlo simulations 3 in the canonical ensemble.These problems in achieving ergodicity arise in a Monte Carlo simulation when the underlying potential energy surface has disconnected but important regions of configuration space separated by significant barriers.At low temperatures the large barriers can trap the system in one region of space so that other important parts of space are either neglected completely or not sampled with sufficient frequency.In the jwalking method, Monte Carlo moves are attempted to all important regions of configuration space using a series of external ergodic distributions generated at high temperatures.The original j-walking method developed for classical Monte Carlo simulations in the canonical ensemble 1 has been extended to quantum simulations 2 and enhanced using histogram techniques. 4Parallel versions of the j-walking method 5,6 are also available enabling ergodic simulations for increasingly large systems.
In the current work we demonstrate that the j-walking method can be extended to microcanonical simulations, and we show that difficulties can be at least as troublesome in the microcanonical ensemble as in the canonical ensemble.Our approach to microcanonical simulations uses methods developed by Pearson et al. 7 and applied to physical systems [8][9][10] including a recent study of Lennard-Jones clusters. 11The contents of the remainder of this paper are as follows.In Sec.II we derive the key expressions for simulations in the microcanonical ensemble.We then show how j-walking can be applied to insure ergodic simulations using the derived expressions.In Sec.III we provide the computational details, and in Sec.IV we apply the methods to 7-and 13-particle Lennard-Jones clusters.We calculate both the temperature and heat capacity and demonstrate the difficulties in attaining ergodic simulations with ordinary Metropolis Monte Carlo methods.We summarize our findings in Sec.V.

A. Basic formulas
In this section we develop the key expressions used for microcanonical simulations.Although the results of this section are not new, we include the discussion to make the j-walking approach to microcanonical simulations as clear as possible.
We begin by expressing the classical microcanonical density of states for a system having energy E, number of particles N and occupying volume V of coordinate space ⍀͑E,V,N͒ϭ 1 In Eq. ͑1͒ H(p,r) is the classical Hamiltonian, and p and r, respectively, label the momenta and coordinates of the particles in the system.Using Eq. ͑1͒ it is easy to demonstrate the standard Laplace transform connection between the microcanonical density of states and the canonical partition function

͑2͒
where ␤ϭ1/k B T with T the temperature and k B the Boltzmann constant.Using the expression for the classical canonical partition function with U(r) the potential energy, the microcanonical density of states can be obtained from the Mellin inversion integral where ⌰(x) is the step function and ⌫(x) is the gamma function. 13Microcanonical averages for any coordinate dependent variable can then be obtained by averaging with respect to the unit normalized density function For example, the average kinetic energy of a system can be expressed ͗K͘ϭ

͑8͒
where ͗͘ denotes an average with respect to Eq. ͑7͒.
A quantity related to the microcanonical density of states is the phase space volume given by ⌽͑E,V,N͒ϭ 1 Using Laplace transform methods similar to those leading to Eq. ͑5͒, it is easy to show that Alternatively, Eq. ͑11͒ can be verified by using Eq.͑10͒ to produce Eq.͑5͒.The entropy in the microcanonical ensemble can be defined by Sϭk B ln ⌽. ͑12͒ In the thermodynamic limit, the definition given in Eq. ͑12͒ agrees with textbook definitions of the entropy in terms of ⍀(E,V,N) ͑to order 1/N). 7By using Eq.͑12͒ as the defining relation for the entropy, for any N the thermodynamic identity (‫ץ‬S/‫ץ‬E) N,V ϭ1/T leads to the expression ͗K͘.

͑14͒
Consequently, we can identify the average temperature in the microcanonical ensemble with the average kinetic energy for any value of N. By direct differentiation of Eq. ͑13͒, the heat capacity C V ϭ(‫ץ‬E/‫͗ץ‬T͘) N,V can be shown to be given by . ͑15͒ Unlike the manifestly positive expression for the heat capacity in the canonical ensemble where ͗͘ T denotes a canonical average, Eq. ͑15͒ is not guar- anteed to be positive for small N.In the microcanonical ensemble, only in the thermodynamic limit can the heat capacity be expressed in terms of a fluctuation ͑of the kinetic energy͒ 7,14 and necessarily be positive.

B. Microcanonical j-walking
To evaluate coordinate dependent properties in the microcanonical ensemble, averages can be evaluated using the weight function given in Eq. ͑7͒.In the Monte Carlo method 15 such averages are determined by executing a random walk where the points in configuration space are visited with probability P E (r).The walk is generated by performing trial moves with some unit normalized distribution T(r 0 →r n ) where r 0 is the current configuration of the system and r n is some trial configuration.Detailed balance provides a sufficient condition for the walk to visit space ͑asymptoti-cally͒ with probability P E (r).Detailed balance can be satisfied 15 by accepting the trial moves with probability pϭmin ͭ 1, P E ͑r n ͒T͑r n →r 0 ͒ P E ͑r 0 ͒T͑r 0 →r n ͒ ͮ .

͑17͒
In the Metropolis method 3 for multiparticle systems, the moves are attempted for one particle at a time with a uniform distribution.If an attempt is made to change the x-coordinate of particle i from x i,0 to x i,n the Metropolis trial probability is given by For the case of such Metropolis walks, the T-factors cancel in the numerator and denominator of Eq. ͑17͒.
In practice such Metropolis walks can cause ergodicity problems that can be expected to be particularly severe in microcanonical simulations.If two regions of configuration space are separated by energy barriers that are greater than the system energy, a walk begun in one region may never reach other important regions of configuration space.As in canonical simulations a cure for such ergodicity problems is possible by executing a high energy walk known to be ergodic and storing the resulting configurations to an external distribution.In microcanonical simulations, the energy rather than the temperature is used as the parameter to generate information about important regions of the underlying potential energy surface.
In a microcanonical j-walking simulation an initial high energy, E h , is chosen so that a Metropolis walk at that high energy can be expected to be ergodic.During the high energy Metropolis walk, configurations are captured periodically and stored to a file.The configurations so chosen need to be sufficiently separated in the Metropolis walk so that correlations between the configurations are broken.At lower energies, a Metropolis walk is again executed for the majority of Monte Carlo points.Additionally, a fraction P J of the moves are attempted to configurations in the high energy distribution, P E h (r); i.e., we take The configurations are chosen from the high energy distribution at random as a further method to break the correlations between points in a Metropolis walk.Because the high energy distribution can be expected to be ergodic and covers all important regions of configuration space, the jumps to the high energy distribution enable the simulation to overcome any energy barriers that separate the important regions of space.To insure detailed balance attempted moves to the high energy distribution are accepted with probability pϭmin ͭ 1, P E ͑r n ͒P E h ͑r 0 ͒ P E ͑r 0 ͒P E h ͑r n ͒ ͮ .

͑20͒
As in canonical j-walking, the procedure discussed in the previous paragraph is useful provided an adequate fraction of attempted moves to the high energy distribution are accepted.When the energy separation between E and E h becomes large, the fraction of accepted moves can become low.To avoid the low acceptance rate, additional external distributions at increasingly low energies are stored for j-walking to still lower energies.When the storage requirements become prohibitive, parallel implementations of the j-walking algorithm can be used. 5In the microcanonical ensemble the procedures for generating multiple distributions are identical to the canonical procedure, and we provide no additional details here.Fuller accounts can be found in the original canonical literature. 1,2

III. COMPUTATIONAL DETAILS
To investigate the utility of the j-walking scheme for microcanonical simulations, we study the energy as a function of calculated temperature and heat capacity as a function of calculated temperature of Lennard-Jones clusters containing 7 and 13 particles.The potential surfaces for Lennard-Jones clusters have many minima, and it is convenient to refer to these minima as isomers.The isomers are separated by energy barriers that can be significant, and ergodicity problems associated with simulations of clusters have been discussed extensively. 16,17The ergodicity problems are particularly serious in ranges of temperatures or energies associated with phase changes.Phase change regions in clusters occur over a range of temperatures and energies where isomerization transitions become important.The phase change behaviors that have often been associated with cluster melting, have been characterized by many differing properties depending on the ensemble used.In molecular dynamics simulations cluster melting has been discussed in terms of fluctuations in the average kinetic energy between distinct rigid solidlike and fluid liquidlike phases. 18In canonical simulations heat capacity anomalies have often been discussed in terms of cluster melting phenomena. 17Phase change regions in Lennard-Jones clusters have proved to provide a fruitful ground for the development and testing of the canonical j-walking algorithm.We find these Lennard-Jones clusters to be similarly useful in microcanonical jwalking simulations.
We describe the Lennard-Jones clusters with the standard potential model where r i j is the distance between particles i and j, and ⑀ and are the standard Lennard-Jones energy and length parameters, respectively.Because we consider energy ranges above the threshold for evaporation, we also enclose each cluster within a spherical hard-wall constraining potential about its center of mass.We choose the radius of the constraining potential to be 4 for both 7-and 13-particle Lennard-Jones clusters, the same value used by Frantz 17 in studies of Lennard-Jones clusters using canonical j-walking.
In the results that follow we compare the heat capacity of the 7-particle cluster and the energy of the 13-particle cluster as a function of temperature obtained using ordinary Metropolis Monte Carlo methods and j-walking methods.For each case the Metropolis calculations are initiated from the lowest energy isomer ͑the pentagonal bipyramid for the 7-particle cluster and the icosahedron for the 13-particle cluster͒, and each energy point consists of 10 7 warm-up moves followed by 10 7 Monte Carlo moves where data are accumulated for the 7-particle cluster, and 10 8 warm-up moves followed by 10 8 Monte Carlo moves where data are accumulated for the 13-particle cluster.At the end of each Metropolis Monte Carlo calculation, the final configuration is stored in a file and used as the initial configuration for the next higher energy point.The temperature is evaluated using Eq.͑13͒, and the heat capacity is evaluated using Eq.͑15͒.
The j-walking calculations are initiated from a random configuration of atoms confined within the constraining volume.The initial energy for the 7-particle cluster is Ϫ10⑀, and the initial energy for the 13-particle cluster is Ϫ24⑀.At the initial energy a Metropolis walk consisting of 10 8 moves is executed, and an external distribution is generated containing 10 5 elements for j-walking from lower energies.At lower energies the calculations consist of 10 6 Metropolis Monte Carlo moves used to warm-up the system, followed by 10 6 Monte Carlo points with jumps made to the external high energy distribution with a probability P J ϭ1/10.When the acceptance rate for attempted jumps drops below 10%, a new external distribution is generated consisting of 10 5 elements.This external distribution is generated using j-walking to the previously created high energy distribution.This process is continued to the lowest energy of the calculation.The entire j-walking procedure is repeated ten times for the 7-particle cluster and 100 times for the 13-particle cluster so that the total number of j-walking and Metropolis Monte Carlo moves are the same.

IV. RESULTS
The computed heat capacity as a function of the computed temperature for the 7-particle cluster is shown in Fig. 1.In Fig. 1 the left hand panel displays the results of the Metropolis Monte Carlo simulation and the j-walking results are given in the right panel.The error bars represent two standard deviations of the mean.The substantial difference in the errors between the Metropolis and j-walking results is evident.Evident from the small graph inset in the left panel is a systematic difference between the two results at low temperatures ͑less than k B T/⑀ϭ0.1).In the inset, the lower graph does not include j-walking whereas the upper curve does include j-walking.The heat capacity determined from the j-walking calculation is consistently higher than the heat capacity found from the Metropolis calculation except for the lowest calculated energies.By examining the structural features associated with configurations obtained from both calculations, we can identify each configuration with an isomer.The structural comparison algorithm, which has been used recently to facilitate the search for isomers in nickel clusters, 19 has been modified to compare a walking configuration to each known isomer.The frequency of occurrence of the isomer is determined by a criterion that uses the lowest average atom-to-atom distance between the known isomer and the walking configuration.At the end of the walk one obtains the probability for the system to be in the ''vicinity'' of a given isomer.This approach is similar to an inherent structure study 20 where the configurations of a fluid are quenched to the nearest isomer.Using these methods, we find that the Metropolis configurations are always identified with the lowest energy isomer ͑the pentagonal bipyramid͒ for the temperature range graphed in the small box of the left panel of Fig. 1, whereas some of the configurations have the structures of the higher isomers in the j-walking simulation at these low temperatures.The cluster isomerizes at an energy lower than predicted by Metropolis Monte Carlo, and the small error bars computed at low temperatures in the Metropolis calculations are artificially low.When the system begins to isomerize in the Metropolis calculation, the computed errors in the heat capacity become large and the computed random errors in the temperature are significantly larger than the resolution of the graph.In the j-walking calculation the errors in the computed temperatures are smaller than the resolution of the graph.
The qualitative shape of the heat capacity as a function of temperature for the 7-particle cluster does not match well the heat capacity curve obtained by Frantz 17 using canonical j-walking methods.In the microcanonical ensemble the heat capacity has a well defined maximum at k B T approximately equal to 0.15⑀.In the canonical ensemble the heat capacity has a plateau at about k B Tϭ0.2⑀.In both the canonical and microcanonical ensembles, C V /(Nk B ) approaches 2.57 as T→0, the equipartition result.
Figure 2 displays graphs of the energy as a function of the calculated temperature for the 13-particle Lennard-Jones cluster.The left panel represents the results of the Metropolis Monte Carlo simulation and the right panel represents the j-walking results.The statistical error is clearly reduced significantly by the j-walking method.The large slope at k B T/⑀ϭ0.28 is similar to that seen in other microcanonical simulations on this system, 21 and the large slope region has often been associated with cluster melting.As might be expected from the importance of isomerization transitions in FIG. 1.The heat capacity as a function of temperature for the 7-particle Lennard-Jones cluster.The left hand panel represents the Metropolis results and the right hand panel represents the j-walking results.Each calculation includes 10 7 Monte Carlo points, and the significant reduction in the statistical noise by the j-walking method is evident.The small inset in the left panel presents both the Metropolis data-͑lower portion͒ and the j-walking da-ta͑upper portion͒ at the lowest calculated temperatures.The systematic deviation is a result of low energy isomerization transitions that are not captured by the Metropolis calculation.the melting region, the difference in statistical errors between the Metropolis Monte Carlo results and the j-walking results are greater in the melting region than for other regions of temperature.We do not display the heat capacity for this system, because the slope of the caloric curve in the melting region is too large to obtain statistically meaningful results, even for the large number of Monte Carlo points used in the calculation.

V. DISCUSSION
We have demonstrated that j-walking methods can be applied to insure ergodic simulations in the microcanonical ensemble using methods that are similar in spirit to the original j-walking method in the canonical ensemble.The difficulties in achieving ergodicity in microcanonical simulations of Lennard-Jones clusters are pronounced and qualitatively greater than difficulties that have been seen in canonical simulations.In the canonical ensemble at any temperature there can be events that can overcome barriers between separated regions of configuration space.In the microcanonical ensemble regions of space separated by barriers greater than the available energy are completely inaccessible to each other.We speculate that the total lack of connectivity between regions of space in the microcanonical ensemble may make it particularly difficult to achieve ergodicity.
The simulations discussed in this work are true microcanonical simulations in the sense that we have made no effort to impose conservation of the total linear and angular momenta as would be imposed naturally in a molecular dynamics simulation on these systems.The imposition of such conservation constraints represents a straightforward modification of the results presented here.

FIG. 2 .
FIG. 2. The energy as a function of temperature for the 13-particle Lennard-Jones cluster.The left panel represents the Metropolis Monte Carlo results and the right panel represents the j-walking results.Each calculation consists of 10 8 Monte Carlo points.