J.-M Soudan

) M Basire

J.-M Mestdagh

C Angelié

A new Monte Carlo method for getting the density of states of atomic cluster systems

A novel Monte Carlo flat histogram algorithm is proposed to get the classical density of states in terms of the potential energy, g(E p ), for systems with continuous variables such as atomic clusters. It aims at avoiding the long iterative process of the Wang-Landau method and controlling carefully the convergence, but keeping the ability to overcome energy barriers. Our algorithm is based on a preliminary mapping in a series of points (called a σ -mapping), obtained by a two-parameter local probing of g(E p ), and it converges in only two subsequent reweighting iterations on large intervals. The method is illustrated on the model system of a 432 atom cluster bound by a Rydberg type potential. Convergence properties are first examined in detail, particularly in the phase transition zone. We get g(E p ) varying by a factor 10 3700 over the energy range [0.01 < E p < 6000 eV], covered by only eight overlapping intervals. Canonical quantities are derived, such as the internal energy U(T) and the heat capacity C V (T). This reveals the solid to liquid phase transition, lying in our conditions at the triple point. This phase transition is further studied by computing a Lindemann-Berry index, the atomic cluster density n(r), and the pressure, demonstrating the progressive surface melting at this triple point. Some limited results are also given for 1224 and 4044 atom clusters.

I. INTRODUCTION

The density of states (DOS) g(E) of a system is a key function, allowing to compute all thermodynamic canonical quantities (free energy F, internal energy U, heat capacity C V , etc.), at any temperature T by simple integration over the energy E. For example, the partition function, from which all other quantities can be derived, is (at the inverse temperature β = (k B T) -1 )

Z(β) = +∞ -∞
dEe ln[g(E)]-βE .

(1)

g(E) is also used in the Rice-Ramsperger-Kassel-Marcus statistical theory of reaction rates, in which the microcanonical dissociation rate is given by

K D (E) = G(E -E d ) hg(E) , ( 2
)
h being the Planck constant, G(u) the integrated DOS of the transition state up to the energy u, and E d the dissociation energy of the system. This formula has been applied to the unimolecular dissociation of molecules 1,2 and clusters 3,4 using Monte Carlo (MC) computations of g and G.

Given this key role, various Monte Carlo schemes have been developed in the last 20 years to determine g(E), such as the multicanonical technique, 5 the entropic sampling, 6 the broad histogram method, 7 and the Wang-Landau (WL) sampling. 8,9 A common feature of these algorithms is to use the Metropolis algorithm with a bias b(E), in which moves a) Electronic mail: jean-maik.soudan@cea.fr.

from energy E a to E b are accepted with the probability

P (a → b) = min 1, b(E a ) b(E b ) . ( 3
)
Consequently, in the standard Metropolis algorithm, the bin interval δE is populated according to an histogram H(E), given by

H (E) ∝ g(E) b(E) . ( 4
)
A flat histogram is found when the bias matches the true density of states (b(E) = g(E)

). An advantage of the flat histogram methods is that the MC point circulates in a large energy range, hence overcoming more easily the barriers in complex energy landscapes.

The WL algorithm has received a keen interest because it is self-adaptive. The initial bias g WL (E) is usually chosen constant and is modified simultaneously with the histogram of visits H WL (E) in terms of a modification factor f. H WL (E) converges to a flat histogram through, typically, 20-25 sequences of MC walks, each one corresponding to a value of the parameter f which starts around e 1 and is forced to tend towards 1, reducing consequently the noise. Then, g WL (E) converges towards the exact density of states. This algorithm has largely been applied to systems with discrete variables (spins), but also, to a lesser extent, to systems with continuous variables, such as Lennard-Jones (LJ) clusters of various sizes: LJ 250 (Ref. 10) (for which g(E) was varied by a factor 10 500 ), LJ 31 , 11 LJ 38 , 12 and LJ 216 . 13 The classical density of states is usually considered, but the quantum one has also been tackled with the WL algorithm. 14,15 These previous works about clusters refer to nanothermodynamics, a fascinating field in which surface effects and non-extensivity play an important role. [16][17][18][19] Particularly, transition metal clusters are crucial ingredients for the catalytic growth of carbon nanotubes, 20 and their thermodynamic state is hardly accessible to experiment and controversial. To tackle theoretically such systems with continuous variables of a few hundreds of atoms or more, tightened by a realistic potential, is an important but very challenging program. The sequences of the WL method slow down when the parameter f approaches 1, and it happens that the convergence is not well controlled, particularly in a phase transition zone. Such convergence problems have been mentioned for LJ clusters. 11,12 In Ref. 10, the WL sequence has been ended by a standard Metropolis walk in view of correcting the density of states obtained. It is also worth mentioning that Kim et al. 21 proposed an adaptive mapping in terms of a temperature T(E) to improve the convergence of the WL f-cascade.

The present work provides a novel approach to derive g(E), particularly suited to systems with many degrees of freedom and with interacting potentials more involved than the LJ modelization. This new method converges in a well controlled way, still acting within the framework of MC flat histogram methods. Its originality is to get quickly a very good initial guess of the bias b(E) through a new type of mapping function σ (E p ) obtained at a few tens discrete points of the potential energy. As a result, b(E) converges in only two reweighting iterations, by MC walks on a few overlapping ranges. This study focuses on the methodology and convergence of the algorithm. This is why we choose as an illustration the model system of a 432 atom cluster of 2 nm diameter, the interaction being given by the simple 6 parameter Rydberg type pair potential, for which billions of MC steps can be performed in one day. Some limited results are also given for clusters of 1224 and 4044 atoms.

This paper is organized as follows. In Sec. II, we explain the basis of our algorithm getting the density of states g(E p ) of the potential energy in the configurational space, in terms of the mapping function σ (E p ). In Sec. III, we illustrate briefly the method on a Debye (i.e., harmonic) cluster, and in detail on the Rydberg cluster, discussing particularly the convergence qualities. In Sec. IV, the canonical quantities U(T) and C V (T) are first derived by integration, to put in evidence a phase transition. This one is further analyzed using the Lindemann-Berry index, the density n(r) of the cluster, and finally a computation of the pressure as a function of the temperature. These analyses demonstrate clearly that this phase transition is a solid-liquid one in the presence of a gas, i.e., at the triple point. In Sec. V, the variation of the mapping function σ (E p ) is explained, and a comparison of the performances with the WL algorithm is made. The conclusion of Sec. VI paves the road for treating a system with more realistic interactions than the Rydberg potential, such as the modified embedded atom method (MEAM) potential. 22,23

II. METHOD TO GET THE DENSITY OF STATES

In this section, we consider the configurational space of the system in which the potential energy E p ( -→ r 1 , . . . , -→ r N ) is computed in terms of the positions -→ r i of the N particles. The Metropolis algorithm is used with standard translational moves in this space, and a fixed bias b(E p ) determined by the method explained below. It provides the density of states g(E p ) associated with the potential energy, whereas g k (E k ) for the kinetic energy is known analytically. Given the separability of the total energy as E = E k + E p , canonical quantities at an inverse temperature β can be deduced by independent integrations over E k and E p , respectively, thermodynamic functions being additive at a given β (see Sec. IV). Thus, only computations relative to E p and corresponding values (g(E p ), σ (E p ), . . .) need to be conducted and are described in the following.

A. Relevance of σ (E p ) to express the DOS g(E p )

This section focuses on quantities ln [g(E p )] (DOS) and ln [b(E p )] (bias), which correspond to the configurational entropy of the system S(E p )/k B . 7,18,19,21 They can be expressed as

ln[g(E p )] -ln[g(E ref )] = E p E ref duβ(u), (5a) ln[b(E p )] -ln[b(E ref )] = E p E ref duβ b (u), (5b)
where the function

β(E p ) (respectively β b (E p )) is the deriva- tive of ln [g(E p )] (respectively ln [b(E p )]
) and appears as a configurational microcanonical inverse temperature. E ref is an arbitrary reference energy. We found convenient and illuminating to express this function β (respectively β b ) under the form

β(E p ) = σ (E p ) E p -E 0 z , ( 5c
)
where E 0 z is in principle the absolute minimum of the potential energy, and σ (E p ) (respectively σ b (E p )) is a mapping function, which is an approximately extensive quantity with no dimension. In fact, Eq. (5c) can be used with an origin E z different from the absolute minimum E 0 z . It must be emphasized that the density of states g(E p ) and the function β(E p ) are quantities independent of the origin E z , whereas σ (E p ) depends on it. This dependence and the importance to determine

E 0 z is discussed hereafter in Sec. II B 1.
The justification of the mapping of Eq. (5c) results from the following model. For a separable potential with N d degrees of freedom of the form

E p = N d i=1 |x i | n i , ( 6a
)
the mapping function is expressed as

σ = N d i=1 1 n i -1, (6b)
and is therefore strictly constant. When the system is in the configurational vicinity of the absolute minimum of the potential energy basin, an harmonic expansion is valid and corresponds to the particular case n i = 2. Then, when E p tends towards this bottom E 0 z ,

σ (E p ) → N d 2 -1. (7)
Note that the application of Eq. ( 7) for a cluster of 432 atoms with 6 frozen degrees of freedom gives N d = 1290 and σ (0) = 644. This is the typical magnitude of σ for the system studied in Sec. III. The previous separable model explains why the mapping function σ (E p ) is more relevant than the parameter β(E p ): (i) Equation (7) shows at low energies that σ (E p ) is finite, provided E 0 z has been determined, while β(E p ) is diverging to +∞. (ii) At higher energies, σ (E p ) has almost everywhere a rather slow variation in terms of E p (except in a phase transition zone). Equations (6a) and (6b) help to understand the variation of σ (E p ) using a modelization of the mean potential perceived by each atom at the energy E p . This discussion is further pursued in Sec. V A. (iii) The two previous properties allow to get a very good representation of σ (E p ) by interpolation (or even extrapolation) between a few points (see first step of the algorithm). (iv) In a phase transition zone, σ (E p ) exhibits a jump on a narrow energy range and displays, therefore, strikingly this phenomenon, while β(E p ) shows only the so-called "back-bending" (an oscillating variation). This is analyzed in Sec. IV A. (v) As discussed in the preliminary step of the algorithm, far from being a limitation of our method, the determination of E 0 z , by the dependency it introduces in σ (E p ), is a powerful self-consistent way of checking, through Eq. ( 7), that no lower basin is reached during all the MC simulations.

It is to stress that all the previous equations and properties are equally valid for b and σ b .

B. Algorithm

The method used to find the mapping function σ (E p ) is implemented in several steps, which are described in detail in the following. They can be summarized as follows:

1. Preliminary step: determination of the minimum configurational energy E 0 z . 2. First step: use of a local linear bias in a series of narrow potential energy ranges to get a first guess σ (1) b (E p ) of the mapping function over an energy grid. Each value is deduced from the Gaussian form of the accumulating energy histogram and interpolated over the energy range to give σ (1) b (E p ). The corresponding bias function b (1) (E p ) is computed from Eqs. (5b) and (5c). 3. Second step: standard Metropolis exploration on large intervals in potential energy, with b (1) (E p ) as an input for the Metropolis rule. A corrected mapping function σ (2) b (E p ) is derived from the histograms H (1) and the mapping function σ (1) b (E p ). Then, a new bias b (2) (E p ) is deduced. 4. To increase convergence and ergodicity, it can be necessary to iterate this last step using b (2) (E p ) in a new Metropolis exploration. Overlapping intervals can now be extended. b (3) (E p ) thus obtained converges now to the real g(E p ).

Preliminary step

The Wales technique 24 is used to search the absolute minimum of E p ( -→ r 1 , . . . , -→ r N ). Obviously, it is not possible to guaranty that the true absolute minimum has been found, whatever the duration of the search. It is, thus, possible that during a longer MC trajectory, a deeper basin could be found and explored. This situation was in fact met in our first simulations: after a too short minimum search, we observed that, during the MC walk, the mapping function σ (E p ) exhibited a sudden fall (see the first step of the method). This provides a very clear signature of a transition to a new and deeper basin. This fact can be explained by the following expression relating two different energy origins E z :

β(E p ) = σ z0 (E p ) E p -E 0 z = σ z1 (E p ) E p -E 1 z . ( 8
)
β(E p ) (Eqs. (5a) and (5c)) does not depend on the choice of the energy origin. The minimum found in the exploration, E 1 z , is a priori larger than the absolute minimum E 0 z . The two functions σ z0 (E p ) and σ z1 (E p ) corresponding to these values are related by Eq. ( 8). Then, if the MC point has explored ergodically, the true minimum basin, the left fraction expressing β(E p ) tends to a finite value when E p tends to E 1 z , and consequently, σ z1 (E p ) must tend to zero in the right fraction, so that the property of Eq. ( 7) is strongly violated. In this case, a new minimum search can be performed from a point extracted from the new basin explored during the MC walk. Hence, the value of the minimum is changed. If the sudden fall of σ (E p ) mentioned above does not occur anymore during the remaining MC explorations, this is the proof that no lower basin has been reached. It is interesting to note that Eqs. (7) and (8) can be used to deduce E 0 z with a very good accuracy from the bad value E 1 z . Choosing a value E p a little above E 1 z and enforcing σ z0 (E p ) to check Eq. ( 7) allows to estimate E 0 z knowing σ z1 (E p ).

It must be stressed again that our algorithm could be implemented without this preliminary search of the absolute minimum, but in this case the correct energy range of the potential energy would not be identified and the property of slow variation of σ (E p ) would be lost.

First step: Linear bias

This step is the key trick of the method. The principle is to map the σ (E p ) function at a series of points E p i , using a series of Metropolis explorations with a mapping function σ L b (E p ), linear in terms of energy, from which a bias b L is derived. Each exploration spans an interval beginning at E α and of half-width α , the two-parameter linear bias being

σ L b (E p ) = σ L b (E α ) + B α • (E p -E α ). ( 9
)
Within the hypothesis of an ergodic exploration, the corresponding histograms H(E p ) is accumulated in this interval according to Eq. ( 4). H(E p ) is, therefore, extremum at the value E p i where the derivative of ln [g/b L ] is zero. Then, the following equality is satisfied:

σ (E p i ) = σ L b (E p i ). (10)
A second-order expansion of ln [g/b L ] around this point gives the expression

H (E p ) ∝ exp[-A(E p -E p i ) 2 ], (11a)
with

A = B α -σ (E p i ) 2E p i , ( 11b
)
determining the width of the function (σ is the derivative of σ with respect to the potential energy, and the energy minimum has been fixed to zero). When the condition If the MC exploration were fully ergodic at this step, the determination of the mapping function σ (E p ) would be achieved. In fact, we emphasize that the MC trajectory is limited to a narrow range of energy, and the system has a limited ability to hop from one energy basin to another. Then, full ergodicity is not reached and Eq. (10) gives only a first (but in fact very good) approximation σ (1) b (E p ), which is used to compute the bias b (1) (E p ) introduced in the next step. The remarkable accuracy of this first step, in spite of the use of narrow energy ranges, is the key for the efficiency of our algorithm.

B α > σ (E p i ) (11c) is satisfied, E pi is a maximum (= E max )
Note that choosing the particular value σ L b (0) = 0 leads to a canonical bias, with a slope corresponding to the inverse temperature, defined at the point E pi by

B T (i) = (k B T i ) -1 = σ L b (E p i ) E p i . ( 12
)
The present two-parameter bias allows a much greater flexibility than the canonical one. This is necessary to get a good accuracy on the determination of σ (1) (E p i ) over all the energy scale, and to map correctly the phase transition zone. The conditions that the parameter B α (and correlatively σ L b (E α )) must satisfy to optimize the accuracy are detailed in Appendix A. The latter sums up the main conclusions. The first condition that B α must meet is Eq. (11c) which provides a lower limit to it. It is constraining only in the phase transition zone where σ is high. Then, the half-width α must be chosen to be a fraction of the energy at the middle of the interval to monitor a sufficiently well-localized zone. Moreover, the Gaussian must be well drawn, centered around the center of this interval. This gives for B α a lower limit much smaller than B T (Eq. (A2c)), which leaves a large possibility of choice. Finally, the determination of σ p using Eq. ( 10) is made by calculating the energy E max where the histogram is maximum. This is performed by fitting the region of this maximum by a third-order polynomial decomposition of σ L b (E p ). A third-order expansion around E max is necessary to take into account a slight asymmetry of the histogram, compared to a Gaussian. This asymmetry contribution can modify the value of σ by several units in some zones and cannot be neglected (Eqs. (A5b) and (A5c)). We have also performed an analysis of the MC fluctuations giving the statistical accuracy δσ (Eq. (A10)). It is interesting to stress that what is important for the flatness of the histograms obtained at the following step is the variation of σ over an energy range and not its absolute value. This is analyzed in Sec. II B 3.

Second step: Metropolis exploration

At this step, a standard Metropolis exploration is done with the previous bias b (1) (E p ). In contrast to the previous step, it is now performed on a series of large potential energy overlapping ranges

[E inf α , E sup α ]. These intervals, indexed by α, have an expansion R α =E sup α /E inf α .
Hence, we now take benefit of the circulation of the MC point in these large intervals to overcome the energy barriers and approach ergodicity. The extension of each range is nevertheless first limited by the need to make a sufficient number of round trips between the extremities of the interval, in view of obtaining good statistics. This restriction is related to the need of a sufficiently high Metropolis acceptation rate, which limits the translational moves. It is stressed that this limitation concerns any flat histogram method using independent intervals. Indeed, only replica exchange methods between several intervals, or concerted movements speeding the diffusion on the energy scale, can extend the energy spanned by the MC point, thus allowing double funnel structures of the energy landscape to be overcome (see Conclusion). It must be noticed that the different intervals overlap, so that they can give slightly different mapping functions σ (2) b (E p ) on the common part of two overlapping energy ranges. This is the consequence of the statistical fluctuations δσ and of a non-ergodic exploration in, at least, one interval. In this case, an interpolation was made between the two mapping functions from one extremity of the common range to the other.

In the frame of the ergodic hypothesis, the (in principle) exact value of σ is deduced from the histogram H (1) and from the mapping function σ (1) b (E p ) by a reweighting derived from Eq. ( 4)

σ = σ (1) b + E p • d ln H (1) dE p σ (1) H +σ cd , ( 13
)
σ cd being a correction due to the discretization of energy into bins δE, given in Appendix B. The parameter δE is chosen in such a way that σ cd is generally negligible, which is the case if δE 2E/σ (Eq. (B6)). It is important to note that a systematic discrepancy σ = σ (1) H over an interval of expansion R α gives a relative variation,

H sup /H inf R σ α . ( 14
)
Hence, the mapping σ (1) b (E p ) obtained from the first step is useful, only if the error σ does not exceed a few units. This is actually the case, as shown in the following. The minimization of the H variation over the interval is, however, a second limitation for its extension R α . The application of Eq. ( 13) gives a new and improved approximation σ (2) b (E p ) of σ (E p ), and then a new bias b (2) (E p ). At this step, it is not known yet if a perfect ergodicity has been reached, so if σ (2) b (E p ) has fully converged towards σ (E p ).

Then, a second iteration can be performed, identical to what has just been explained, but with the new estimate b (2) of the bias. With this new bias, the ranges R α can be extended to improve the ergodicity and diminish the number of interval junctions. The extension of the range remains nevertheless limited by the need of a sufficient number of round trips. This last iteration was found to be necessary at low energies, in the solid zone, where the ergodicity was difficult to obtain. On the contrary, in the fluid zone this iteration serves only to check that b (2) gives a flat histogram.

Final remarks

We conclude this section remarking that the previous σ mapping is suited to a system with continuous variables, such as the position coordinates described here in a classical picture. Its adaptation to, for example, quantified spin systems remains to be explored but would need another type of slowly varying parameter.

It is also to notice that the linear bias σ L b (E p ) is akin to the temperature-linear bias used in a recent work, 25 in the context of the replica exchange method. In this work a temperature T α (E) depending linearly on E is used in a set of intervals α to compute a bias w α (E). For a sufficiently high slope, a localized histogram with a Gaussian shape is built. The set of noncanonical ensembles α is used in a replica exchange method. At the end, the entropy S(E) is calculated via the weighted histogram analysis method. Contrary to that, our σ -linear mapping provides directly a first estimate of ln [g(E)] (or S(E)).

III. ILLUSTRATION OF THE METHOD AND DISCUSSION OF THE σ (E p ) CONVERGENCE

A. Computational details

The calculations were performed with a PC Dell T5400 with eight processors of 2.86 GHz. The initial crystal is cut spherically with a radius of 10 Å, in bulk bcc at T = 0 K, corresponding to 432 atoms. The crystal is placed in a spherical cavity of radius R W = 14 Å, leaving 4 Å free around the crystal. An atom is fixed at the center to prevent translation and three other coordinates are frozen to prevent rotation. The same process is also used to cut spherical clusters of 1224 and 4044 atoms, with a single R W = 25 Å sufficiently large to contain either of them.

B. A Debye crystal

As a preliminary test of the method, a Debye crystal has been studied, i.e., atoms bound to the first neighbours by an harmonic potential with a unique constant K ij = 1 eV/Å 2 . In the range [50-100 eV], σ (E p ) = 644.10 has been obtained in 10 9 steps, in very good agreement with the analytical value for an harmonic cluster having N d = 1290 degrees of freedom: σ = 644 (Eq. ( 7)). However, in the range [4000-8000 eV], when atoms collide with the wall, σ (E p ) decreased to 641.8. This is not surprising given the reduction of the accessible configurational space due to the wall (see Sec. V A).

C. A Rydberg crystal

The atomic interaction is described hereafter by a pair potential of the Rydberg type, analogous to the Rose function of the MEAM potential 22,23 =

-E c [1 + a * + d(a * ) 3 ]e -a * , ( 15a
)
a * = α(r/r e -1). (15b) A cutoff function f c (x) = (1 -(1 -x) 2 ) 2 is used for 0 < x < 1, with x = (r c -r)/δr c .
The six parameters are chosen to reproduce approximately the bulk data of iron around the equilibrium bcc structure. Nevertheless, the purpose of this work is not a quantitative study of this system, which cannot be well represented by a pair potential. 26 The parameters are: E c = 0.6371 eV, r e = 2.64922 Å, α = 5.170477, d = 0.05, r c = 4.0 Å, and δr c = 1.144422 Å. With this potential, 2.4 × 10 9 MC steps can be performed for 432 atoms in one day.

The preliminary step of the method presented in Sec. II gave, after a long Wales search (about three days of computing time), an optimized minimum E 0 z = -1434.8245 eV, with 13 other minima in a range of 0.6 eV. Note that an accuracy of 10 -4 eV for E 0 z is necessary, when exploring intervals beginning at 0.01 eV of the absolute minimum, in view of checking Eq. (7). In all the following MC explorations, we never encountered the fall of σ at low energies that can happen after a too short preliminary search, as it was mentioned in Sec. II. This provides high confidence that the absolute minimum has been found. In the following, if not otherwise stated, the origin of the energy is chosen at this absolute minimum (instead of being at the infinite separation of atoms) and, therefore,

E 0 z = 0.
Step 1 is illustrated in Fig. 1, which displays a Gaussian histogram obtained with a σ -linear bias near 350 eV. The full mapping σ (E p ) is obtained with 25 points distributed from 1 to 6000 eV (not shown here). About 5 × 10 9 MC steps were necessary for a good convergence below 200 eV, but above this value, 5 × 10 8 were sufficient. It is important to stress that this preliminary mapping with 25 points has a computational cost corresponding to at most two intervals of step 2, so that it supplies a very good first guess of the σ curve at a remarkably low cost. In particular, it localizes strikingly and accurately the phase transition zone, which will be discussed in Sec. IV. Indeed, this curve shows a large slope of σ (E p ) near 185 eV, the corresponding temperature being quantitatively determined at this stage within 1%.

The value of the energy at the maximum of the histogram presents fluctuations which determine the accuracy of the σ value. We find δσ 1 between 100 and 1000 eV, i.e., typically one-thousandth of the value of σ . This uncertainty is linked to the MC statistics and must not be mistaken with the error relative to the exact σ obtained from an actual ergodic exploration almost achieved in the following steps. δσ is analyzed in Appendix A (Eq. (A10)) and its expression explains correctly the observed fluctuations.

Usually, the starting MC point of each Gaussian is prepared in the following way. The system evolves freely from a point obtained at lower energy towards a higher energy range of interest. The dynamical evolution of the Gaussian in this new interval can be studied, and it is observed that the energy of the maximum increases up to an asymptotical value. Moreover, a striking phenomenon appears in the phase transition zone just below the high slope region and is observed directly on the histogram dynamics. A single Gaussian is first formed as usual, but suddenly a new peak appears at slightly higher energy and grows to become a new Gaussian. This corresponds to a σ jump of several tens. This will be analyzed in more detail in Sec. IV B and interpreted as a transition of the MC point from a solid basin to a liquid one.

Step 2 is performed on the full energy range [0.01 -6000 eV], covered by 12 overlapping intervals. The result is illustrated in Fig. 2(a) for a large energy range [230 -1530 eV], typical of the fluid zone (line "first iteration"). The guess b (1) (E p ) was obtained in this range by 12 points coming from σ -linear bias explorations. The resulting histogram H (1) (E p ) after 50 × 10 9 steps varies between 4 × 10 -2 and 1, implying an average correction σ over the interval by +2.5, i.e., only a few thousandths of the σ mapping estimate. This correction is due to the lack of ergodicity of the previous step. With Eq. ( 13), a new estimate b (2) (E p ) of the density of states is computed.

A second iteration is performed with the bias b (2) (E p ) previously determined, using only eight overlapping intervals. The result is shown in Fig. 2(a) (line "second iteration"), demonstrating that the convergence to a flat histogram is achieved in this interval. This proves the necessity of large energy ranges to explore correctly the configurational space. This fact is also supported by the following observation: when performing step 2 just below the phase transition at 185 eV, and extending the energy range towards low values (50 and next 25 eV), the new correction σ (2) H reaches several units below 100 eV (+7 near 50 eV) at the end of this step. It is the evidence of a lack of ergodicity, which is not surprising for ranges remaining entirely in the solid zone. The circulation of the MC point near the fluid zone is necessary to approach ergodicity and to overcome energy barriers, as it was also shown in Refs. 11 and 12 in the context of the WL method.

Interestingly, the average number of steps to make a round trip between the extremities of an interval α of typically R α 4 is only 1 × 10 7 in the [14-55 eV] range. It raises up to 5 × 10 8 in the [50-160 eV] range just below the phase transition and to 3 × 10 8 in the [230-1530 eV] range just above. It falls down to 7.2 × 10 5 in the high energy range [1400-6000 eV], indicating the easiness of convergence when the phase is very fluid. We have also observed that the average number of MC steps necessary to make a single round trip varies approximately as the square of the interval extension in a given zone (see Appendix A, Eq. (A8)). This is the indication of a diffusion process over the energy scale. 27,28 As shown in Appendix A, the uncertainty of the histograms is inversely proportional to the square of the number of round trips (Eq. (A8)). We conclude from these statements that it is necessary to achieve a compromise: sufficiently large intervals must be defined to approach ergodicity in the solid zone, but not too large to be able to make a few tens round trips, and thus to get a good statistical accuracy of about 10% on g(E p ).

Another striking phenomenon is seen in the narrow [150-220 eV] range around the phase transition at 185 eV. This phenomenon is illustrated in Fig. 2(b) by successive snapshots of the histogram in this range. While the system wanders rapidly inside each half-range below or above 185 eV, 8 × 10 9 steps at least are necessary to switch from one half-range to the other, so that only three round trips are achieved during our 50 × 10 9 steps MC walk, although a flat histogram method is used. This is the evidence of a strong kinetic barrier between the two basins corresponding to the Figure 4 shows the variation of Log 10 [g(E p )] indicating that g(E p ) spans 3700 orders of magnitude when the energy varies over the range [0.01-6000 eV]. Comparing with Fig. 3, it can be noticed that the phase transition is undetectable on the curve Log 10 [g(E p )].

Finally, the curve σ (E p ) was computed for the 1224 and 4044 atom clusters, from the minimum energy up to the liquid zone, in view of obtaining the melting temperature in terms of the cluster size. We also found for these two clusters that it was not possible to make a single round trip between the solid and the liquid zone for an interval overlapping the phase tran-FIG. 4. Density of states g(E p ) plotted for a 432 atom system from 0.01 to 6000 eV in Log 10 scale. sition region. This was the case even for the 1224 atom cluster in 30 × 10 9 steps, and a fortiori for 4044 atoms. Moreover, we observed that a MC configuration placed near the phase transition but in the solid zone was quickly attracted to the liquid zone with no return. This is the sign that the kinetic barrier from the liquid to the solid phase increases with the number of atoms. However, when the interval was completely inside the solid zone or the fluid zone, the MC point could make many round trips between the extremities, whatever the number of atoms in the cluster. It is stressed that, in spite of the impossibility to transit between the solid and the liquid zones, it was possible to accurately determine σ (E p ) for such large clusters. This was achieved by exploring an interval ending just below the phase transition zone (within a few eV), and an interval starting just above. That way, a spline interpolation in the few eV range of this phase transition zone leaves only a small uncertainty.

The resulting σ (E p ) curves are represented in Fig. 5, in units per atom, and show the same general shape as for the 432 atom cluster. The most obvious difference is that the energy of the triple intersection is shifted to higher values with the increase of the number of atoms.

IV. CANONICAL QUANTITIES AND PHASE TRANSITION CHARACTERIZATION

A. Caloric curve and heat capacity

Taking kinetic energy into account, canonical quantities, such as the direct caloric curve, can be computed as

U(β) = E k + E p , with E k = N d /2β, U p = E p = 1 Z p +∞ -∞ dE p E p e ln[g(E p )]-βE p , ( 16
)
Z p being the configurational partition function (see Eq. ( 1)), and β the constant defining the inverse temperature (which must not be confused with the function β(E p )).

The previous equation allows a discussion of the curves shape reported in Figs. 3 and 16) has one maximum at this point. This situation occurs when σ < β, for which Eq. ( 16) is very well approximated by a Gaussian expansion analogous to that of Eq. (11b). It results that U p E β p and the configurational heat capacity C pV (T) = (∂U p /∂T) V is simply related to σ by

σ C pV k B 1 - σ β . ( 17
)
σ follows consequently the variation of C pV since the bracket is positive. When the equation β = β p (E p ) has a triple intersection, the integrant has two maxima and one intermediate minimum. [10][11][12][17][18][19] This situation occurs when σ > β for the central intersection, and reveals a phase transition at the equality of these two maxima. The phase transition zone corresponds, therefore, to a high slope of σ . Then, U(β) jumps from the first maximum to the second, and (∂U/∂β) V diverges to high values, implying the presence of a peak in the heat capacity. It results that σ does not follow anymore the variation of the heat capacity in the phase transition zone, but is rather related to the jump of U. To show evidence of the phase transition zone, the heat capacity is computed in the following by the Einstein relation on the variances:

C V (T ) = (∂U/∂T ) V (18a) = k B β 2 E k 2 + E p 2 . ( 18b
)
E p and E 2 p are computed by numerical integration over E p , knowing g(E p ). The kinetic part is

N d /2k B .
The resulting values for U(T) and C V (T) are depicted in Fig. 6. In the direct caloric curve U(T), a large slope is observed around the phase transition. The leap of 52 eV is associated with the latent heat of liquefaction. The phase transition is more clearly localized on the curve C V (T), with the peak centered at T ph = 2275 K (β = 5.10 eV -1 ), having a FWHM of 56 K.

For the clusters with 1224 and 4044 atoms, the phase transition is found at, respectively. T ph = 2805 and 3271 K. The temperatures T ph of the three clusters are reported in Fig. 7 in terms of N -1/3 at , which is proportional to the inverse of the cluster radius. The three points form a straight line, as expected from previous studies, [29][30][31][32] and the bulk value can be extrapolated to T ph = 4150 ± 50 K. This high value is characteristic of the Rydberg pair potential which is not a realistic representation of metallic compounds because of the lack of an environment dependence 26 (for comparison, T ph = 1811 K experimentally for bulk iron).

B. Lindemann-Berry analysis

Phase transitions can be probed at the atomic scale using the Lindemann-Berry order parameter, 33,34 called L B . Each MC computation is localized in a narrow energy range E using the σ -linear bias at the step 1 of the method. A spatially restricted version L B (R 1 , R 2 ) is used here, to monitor the spherical shell Sh between the radius R 1 and R 2 , and is defined by

L B (Sh, E) = 2 N pair i,j N MC ij r 2 ij r ij 2 -1 1/2 , ( 19a
)
where N MC ij is the number of MC steps for which both atoms ij of a pair are present in the shell Sh. N pair is the sum over all pairs ij of N MC ij . The means are summed over the MC steps α

r 2 ij = 1 N MC ij α r (α) ij 2 , ( 19b
)
r ij = 1 N MC ij α r (α) ij . ( 19c
)
It is important to notice that the L B coefficient is a dynamical one, either in MC or MD computations, and it depends on the number of steps performed. For an infinite number of steps, the full configurational space would be explored. L B would necessarily characterize a diffusion of the pairs in the full accessible volume, with the limit value 0.375 whatever the potential energy. On the contrary, for a finite number of steps, a low value of L B is expected in the solid phase where the movements remain local. Then, the L B coefficient calculated for the core (0 to 9 Å) and the surface shell (9 to 12 Å), after either 2 × 10 9 or 9 × 10 9 MC steps, are shown in Fig. 8. An increase of the number of MC steps simply translates the curves to lower energies of typically 20 eV, the diffusion being more extensive. Below 80 eV, the L B coefficient is characteristic of local moves in a solid phase ( 0.05), both for the core and the surface atoms. At high energies (>200 eV), L B converges to the value 0.375, which characterizes the full diffusion of a fluid phase in all the accessible volume. Between these energy ranges, the L B curves indicate markedly that the surface shell melts more easily than the core. Indeed, the surface melting extends from 80 to 200 eV, while the core melting begins later, near 115 eV. The energy range between 80 and 200 eV, where this behavior is observed, matches the phase transition zone revealed above by the functions σ (E p ), U(T), or C V (T). In fact, the point where the σ curve increases steeply lies near the end of this range, at 185 eV. Hence, the present Lindemann-Berry analysis identifies the solid-liquid transition. Nevertheless, as the wall is located about 4 Å far from the T = 0 K cluster surface, the presence of a gas around the cluster is expected at this phase transition. Hence, three phases would be coexisting, the present calculation rather corresponding to a triple point. Sec. IV C states this hypothesis.

C. Volumic density analysis

In view of clarifying the presence of gas around the condensed cluster, the atomic density n(R 1 , R 2 ) was also monitored in successive shells between the radius R 1 and R 2 , to get a view of the accessible volume. The resulting density n(r) is given in Fig. 9 at different energies. As for the Lindemann-Berry analysis the calculation is made at step 1 of the method in a narrow energy range, defined by a σ -linear bias.

In the solid phase (80 and 135 eV), or in the liquid phase (E p ≥ 235 eV), the presence of a gas of increasing density (in terms of energy) can be observed between 12 and 14 Å. The surface defining the condensed cluster extends from 8 to 12 Å. Then, at energies corresponding to the solid-liquid phase transition zone, from 135 to 235 eV, a gas is actually present around the condensed cluster. The average number of atoms between 12 and 14 Å increases with energy, from 2.5 atoms (135 eV) to 24 atoms (235 eV). Hence, the solid-liquid phase transition discussed in Sec. IV B occurs at the triple point.

It is interesting to discuss the role of the wall radius choice, R W , on this triple point. When R W is smaller than a value R T R1 W 12 Å, the cluster is expanded up to the wall when energy increases, expelling the gas phase before liquefaction can occur. Then, the solid phase is compressed without presence of a gas phase and the triple point is skirted round to upper pressures. In the opposite case where R W is greater than a value R T R2 W , the cluster vaporizes in the initial void next to the condensed phase, when energy increases, before liquefaction can occur, and the triple point is skirted round to lower pressures. The condition of the previous computations corresponds clearly to the intermediate case

R T R1 W < R W < R T R2 W
where the triple point is reached. The exact value of R W and the total number of atoms determine the size of the cluster at this triple point and its precise characteristics.

A striking fact appears in Fig. 9. The position of the density fall, defining the surface of the condensed cluster, plotted in terms of energy, is a fixed point at 10.5 Å up to E p 1100 eV, where the density becomes fully homogeneous inside the accessible sphere. This property can be interpreted by the fact that the liquid-gas equilibrium is kept until the critical point in pressure is reached at energies near 1100 eV. This situation can be explained following a reasoning similar to that developed above for the triple point. Indeed, when R W is smaller than a given value R CR W , necessarily greater than R T R1 W , the dilatation of the condensed phase is dominant over vaporization and the liquid phase reaches the wall, eliminating the gas phase. In the opposite case, R W > R CR W , the vaporization is complete and the cluster surface moves towards the center of the sphere (r = 0), R CR W being necessarily smaller than R T R2 W . The situation of Fig. 9 corresponds to R W R CR W and happens accidentally because of our combined choice of 432 atoms and of a volume V = 11 494 Å 3 . The latitude on the choice of the volume to remain in this particular situation is related to the fluctuation of the atom number in the external gaseous shell (r ∈ [10. 5-14 Å]). This statistical fluctuation is estimated to be ±16 atoms, corresponding to a latitude of about 4% for the total volume. The precise characteristics of this critical point depend now only on the total number of atoms.

D. Pressure study

Although the condition of a fixed wall, infinitely repulsive, is not well suited for a pressure study, the previous statements that the triple (respectively critical) point has been reached at energies around 185 eV (respectively 1100 eV) suggest to study the pressure curve in terms of temperature. The pressure is calculated for energies E p > 200 eV by using the virial theorem

P V = N at k B T + 1 3 V ir,int , (20a)
with V ir,int being the internal virial term for N at atoms in a volume V, expressed as

V ir,int = i -→ r i . -→ F i . ( 20b
)
The virial is computed by MC simulations with a canonical bias. For E p < 200 eV, the accuracy of this method is completely insufficient, owing to large fluctuations of the virial. Then, the atomic density n W near the wall has to be calculated to estimate the gas pressure using the perfect gas law. The correspondence between the two methods is very good in the range [200-400 eV], and the perfect gas law becomes more accurate at low temperatures (hence energies) for which densities are smaller.

In 10, the resulting pressure curve P(T) is plotted in semi-logarithmic scale, at constant N at and V. In the left part, the solid-gas equilibrium at temperatures T < 2275 K is observed, while for T > 2275 K the liquid-gas equilibrium is obtained. T = 2275 K is the triple point, with a pressure of 145 bars, corresponding to an average of 23 atoms vaporized in the shell between 12 and 14 Å. The right part provides the pressure curve dilated up to 120 000 K, where a homogeneous fluid phase is completely settled with a density n = 3.75 × 10 -2 atoms/Å 3 . Then, the gas obeys a nearly perfect gas law (the virial term being only 10% of the kinetic term). The homogeneity of the density is reached at an energy of about 1100 eV, corresponding to a temperature T c = 25 200 K, and a pressure P c = 156 000 bars. This is the critical point for the 432 atom system as discussed in Sec. IV C. For comparison, bulk iron has its critical temperature in the range 9250-9600 K, and its critical pressure in the range 8250-10 300 bars, 35 but we stress again that the Rydberg pair potential is not suited for a metal. The left part depicts the zone of the triple point from 1000 to 3500 K. The right part is expanded up from 3500 to 120 000 K and shows the critical point.

V. DISCUSSION RELATIVE TO σ (E p ) AND g(E p )

A. Discussion of the σ (E p ) variation

A quantitative discussion of the σ (E p ) curve shape (Fig. 3) is outside the scope of this work, but we can qualitatively sketch the following process. At small energies, the configurational point is in the harmonic well of the minimum energy basin, so that σ is close to 644 (Eq. ( 7)). At higher energies, the negative anharmonicity (sub-harmonic potential, n i < 2) and the appearance of excited basins, corresponding to local rearrangements of surface atoms, produce an increase of σ (E p ). This increase of σ (E p ) can be understood qualitatively assuming that each atom coordinate moves at the energy E p in a mean potential roughly represented by the expression of Eq. (6a). The mean potential becoming sub-harmonic, n i decreases below 2 and Eq. (6b) explains the increase of σ (E p ). The phase transition zone is associated with the opening of a completely new basin corresponding to a global rearrangement, with a higher local energy minimum, but looser (i.e., with smaller frequency vibrations and more negative anharmonicities). This results in a sharp increase of σ (E p ) (cf. Fig. 3(a)). The system can transit between the two types of basins thanks to the strong variation of a collective variable acting as an order parameter.

Then, a maximum is reached (σ 1140 for the 432 atom cluster, cf. Fig. 3(b)) and the subsequent decrease can be explained by the following simple and semi-quantitative model. At high energies, the system behaves as a gas of high volumic density, for which the attractive part of the potential becomes negligible and the repulsive part is dominant. Such a gas can be described again in the approximation of the mean potential perceived by each atom. Within these approximations, the virial term can be expressed through an excluded covolume V c (E p ) perceived by each atom in the presence of all the N at -1 other atoms, which corresponds to a van der Waals analysis. 36 At a qualitative level, Eq. (6a) can be evoked again: a stiff repulsive part corresponds to a strongly increasing n i , much greater than 2, and to the fall of σ (E p ).

At a more quantitative level, the covolume is modeled in the frame of the following argument. The repulsive part of the mean potential is assumed to be proportional to 1/(rr 0 ) 1/q , where r 0 is an atomic hard sphere radius, defining a covolume limit V 0 c reached at infinite energy, and 1/q is an exponent linked to the stiffness of the potential repulsive part. The expression of the excluded covolume is easily deduced, once a reference energy E 0 p (fixed at 1500 eV) has been chosen

V c (E p ) = V 0 c + δV 0 c E 0 p /E p q . ( 21
)
V 0 c + δV 0 c is the covolume at the reference energy E 0 p . The excluded covolume V c (E p ) is a decreasing function in terms of E p . V 0 c , δV 0 c , and the exponent q are three parameters which have to be determined. The integrated DOS, G(E p ), is the total accessible configurational space volume, proportional to (V -V c (E p )) N at . The function σ (E p ) can then be derived in this energy zone from the definition of Eqs. (5a) and (5c)

σ (E p ) = q(N at -1) δV 0 c V -V c (E p ) E 0 p E p q -(q + 1), (22)
where V = 11 494 Å 3 is the total volume of the sphere. V 0 c is estimated to be 0.161 × V by assuming a cubic gathering of all the other N at -1 atoms at a distance defined by a* = -2 (Eq. (15b)), i.e., in the highly repulsive part of the potential. It is now possible to reproduce the curve σ (E p ) between 1500 and 6000 eV by fitting the two coefficients: δV 0 c = 0.217 × V and q = 1.46. Then, the decrease of σ (E p ) in this large energy range can be interpreted as due to the limitation of the accessible volume by the highly repulsive part of the potential, and its energy dependence.

B. Comparison of the σ -linear bias method with the Wang-Landau algorithm

In this section, a comparison of the σ -linear bias method with the WL algorithm is performed, focusing on the analysis of convergence properties, either speed or accuracy. Note that this comparison is relevant only for systems with continuous variables. Three intervals are investigated, ranging from the solid zone ([81, 171 eV], with a ratio R α = 2.11) to the liquid phase ([231, 1530 eV], R α = 6.62), through a phase transition zone ([150, 210 eV], R α = 1.40).

The WL algorithm is applied in the standard way, 8, 9 the modification factor decreasing according to the law f → √ f , from f = 4. The process is stopped when f < 1 × 10 -6 . At a given value f, the MC exploration starts from the value g WL (E) obtained at the end of the previous sequence. The first criterion of stopping is that the minimum value (H min ) of the visiting energy histogram (H WL (E)) becomes greater than (1 -WL ) × H WL (E) , with the common choice 1 -WL = 80%. According to the rule suggested by Zhou and Bhatt, 37 a second stopping criterion has to be simultaneously verified: H min > K/ √ ln f (K = 10 is chosen). This criterion guaranties that the ultimate fluctuations limit of H WL (E), 1/ √ ln f , is smaller than H min . It turns out that this second condition was always checked before the first, and is, therefore, not restrictive. Finally, the number of MC steps, N f , is recorded for each f-sequence, as well as the number of round trips over the interval (ν RT ). At the end of the WL procedure, a simple Metropolis sampling is performed to check the convergence, as made in Ref. 10, the final g WL (E) being used as a bias b(E).

It is observed that the first f-sequences give an estimated density ln[g WL f (E)] strongly biased compared to the final estimate, the discrepancy being progressively reduced when f approaches 1, as usual for the WL methods. 38 The use of a decorrelating parameter S, as proposed by Zhou and Bhatt, 37 allows to strongly diminish this bias, but at the price of an important slowing down of the WL process. For the phase transition range, MC samplings were tested with the decorrelating parameters S = 1, 100, and 1000, the histogram being implemented only after S steps. The initial discrepancy of ln [g WL (E)] to the converged density was strongly diminished with the increasing of S, at the beginning of the WL procedure. Nevertheless, the total number of steps of the f-cascade was heavily lengthened from 2.6 × 10 10 for S = 1, to 9.4 × 10 10 for S = 100, and 3 × 10 11 for S = 1000. The decorrelating parameter is, therefore, practically useless and only the S = 1 value is discussed in the following.

The histogram resulting from the last Metropolis sequence (WL + ), in the phase transition zone, is plotted in Fig. 11, compared to that of our algorithm (identical to curve (3) in Fig. 2(b)). It appears clearly that the last f-sequence of the WL algorithm has not converged properly, the population of the energy range above 185 eV being only 4.1 × 10 -2 (corresponding to a logarithmic variation of 3.2 between the two half-zones).

The WL accuracy predicted by Zhou and Bhatt 37 or Lee et al. 9 on ln[g WL f (E)] is √ ln f , i.e., 10 -3 for the histogram of Fig. 11. It is obvious that the convergence obtained is far from this theoretical limit and it is necessary to explain this bad convergence of the WL algorithm in the phase transition zone. Indeed, the previous estimate of the accuracy is a mathematical limit, reached only at an infinite number of steps and is related to the fluctuations of the continuous self adaptation of the bias. In fact, the practical accuracy in the f-sequences turned out to be also limited by the number of round trips ν RT , as in our method. The importance of the number of round trips has been discussed only by few authors, 27,28 but not in the context of the WL algorithm. This role of ν RT on the accuracy is observed strikingly in the phase transition zone ([150-210 eV]), and deserves now a careful analysis. The result of Fig. 11 can be explained by studying the dynamical evolution of the WL histogram H WL (E) and of ν RT . It is observed that the MC point wanders rapidly inside each half-range below and above the critical value (185 eV), but the crossing of the kinetic barrier needs several billion steps. This is the same phenomenon as observed in step 2 of our method. The number of MC steps necessary to make one round trip, N RT , must therefore be compared to the number of steps of each f-sequence, N f . The evolution of N f in term of ln f during the WL process is reported in Fig. 12 for all studied intervals. N f increases strongly when f decreases towards 1 over all energy ranges, with successive stages, starting from N f 10 8 steps at f = 4. For the phase transition zone, N f increases particularly strongly for the last f-sequence, for which three round trips were recorded. This leads to an average number N RT of 6 × 10 9 steps for a single round trip. For this reason, except for the last f-sequence, a single crossing of the kinetic barrier occurs.

Consequently, in all but the last f-sequence: (i) only onehalf of the full zone is visited in a first part of the WL exploration; (ii) then, the second one is populated after the crossing; (iii) finally, the criteria for stopping are achieved, but no back turn to the first zone has occurred. Hence, no statistical equilibration is settled between the two half-zones. This means that a third necessary condition to stop an f-sequence is needed. The answer would be to enforce a sufficient number of round trips between the extremities of an interval.

The last f-sequence of 17 × 10 9 steps with a limited equilibration of three round trips deserves a particular analysis. A characteristic of the WL algorithm is the following logarith-mic relation between the density of states and the histogram for a bin i:

37 ln g WL i = ln(f ) • H WL i -ln[θ i ], (23)
θ i being the initial bias of the f-sequence. It results that the fluctuations of the histogram are amplified exponentially on g WL (E). An estimate of these fluctuations between the two halfzones, in terms of the number of round trips is outside the scope of this work. However, at the end of the first round trip, it can be estimated that the difference of the accumulated histogram of visits between the two half-zones, δH WL i , is about 2 WL N RT /N b , N b being the number of bins in the interval. From the above equation, this implies a resulting discrepancy on the density of states, for

ν RT = 1 ln[g WL ] = ln(f ) • 2 WL N RT N b . ( 24
)
For the subsequent evolution of the fluctuations, it can be guessed that they diminish in terms of the number of round trips until their ultimate limit 1/ √ ln f is reached. Equation (24) can be applied to the case of the phase transition zone with ln f = 10 -6 , WL =0.2, N b = 1200, and N RT = 6 × 10 9 steps. The predicted fluctuation is ln[g WL ] = 2.0, which must be compared to the observed variation of the logarithm in Fig. 11, i.e., 3.2. Equation ( 7) predicts, therefore, the correct order of magnitude of the fluctuations for a small number of round trips. This result is in contrast with the standard Metropolis approach we applied, for which the linear relation of Eq. ( 4) is available, leading to the following expression: ln[g] ∝ 1/ √ ν RT , which is always smaller than 1 and decreases in terms of the number of round trips. On the contrary, many round trips are recorded in the two other intervals (solid, liquid) allowing a much better convergence. Indeed, the histograms obtained during the last Metropolis stage (WL + ) are more flat. For example, in the fluid zone ([231, 1530 eV]), the histogram WL + oscillates between 0.7 and 1. Nevertheless, it indicates that a last correction has to be made according to Eq. ( 4) to get a flatness similar to that obtained with our method (cf. Fig. 2(a)). Consequently, this last Metropolis exploration is needed, even in these intervals, to improve convergence. It is interesting to apply Eq. ( 24) to the fluid zone with ln f = 10 -6 , N b = 13 000, and N RT = 3 × 10 8 steps, which gives ln[g WL ] = 2.3 × 10 -2 at the end of the first round trip. This is much smaller than the discrepancy of 0.3 which has been observed after H T = 3.5 × 10 9 steps, corresponding to 12 round trips. Equation (24) applied with H T instead of N RT would give the correct order of magnitude ln[g WL ] = 0.276. This is probably the indication that the amplitude of fluctuations does not decrease after a small number of round trips in the WL process. It can be concluded that the convergence of the WL algorithm towards the ultimate limit 1/ √ ln f is very slow. As can be noticed for the three intervals in Fig. 12, the first 15 sequences of the WL f-cascade are fairly fast, and give a rough approximation of ln [g(E)]. The comparison with our method can be detailed for the fluid zone. These first 15 sequences totalize 3.5 × 10 9 steps. Unlike this, the σ -linear bias (step 1) provides directly a better guess with a similar computational effort: 12 Gaussians converge in 5 × 10 8 steps each one. The last seven sequences of the WL f-cascade are much slower (14 × 10 9 steps, i.e., 80% of the calculation time). Only the last sequence of 3.5 × 10 9 steps determine the final accuracy of 30% with its 12 round trips, so that the six preceding sequences waste the computational time. In contrast to that, a single Metropolis sampling with the bias b (1) (step 2) of our method gives, in the same computational time of the last f-sequence, a corrected and converged DOS with an accuracy 8%, although in a more pedestrian way.

VI. CONCLUSION AND PERSPECTIVES

A novel method is described in the present paper, providing the density of states of assemblies of several hundreds to several thousands of atoms. A local probing by a σ -linear bias has been shown to give a very good starting point, at a low computational cost, to obtain at the next steps flat histograms over a small number of large intervals, with only two subsequent iterations. Large clusters up to 4044 atoms, i.e., of a 4.2 nm diameter, have been treated. The present article reports a detailed analysis of the convergence, accuracy, and optimization of the method, using a fairly simple potential in order to be able to perform extensive MC sampling. This analysis is explained in Appendix A and gives informations of general interest for flat histogram methods. Section IV shows how this method, combined with various indicators, provides a powerful tool to analyze the thermodynamic state of a cluster. The curve σ (E p ) obtained in such a way presents a characteristic shape linked to the physical and thermodynamical state of the cluster, particularly in the phase transition zone. It would deserve to be analyzed in the framework of a model valid from the solid to the fluid zone, more sophisticated than that of Sec. V A.

It remains nevertheless a difficulty of convergence below 100 eV, as mentioned at the end of Sec. III, due to an incomplete ergodicity in the solid phase. This has occurred in spite of the exploration by a flat histogram method of large intervals extended up to the phase transition zone at 185 eV. An elegant way of bypassing this problem, at no computer time expenses, would be to combine the present strategy with the parallel tempering/exchange replica approach, 25,39,40 i.e., by exchanging the MC points between overlapping ranges. This combination would bring our method in kinship with the recent work of Kim et al., 25 who have made parallel tempering between Gaussian Metropolis explorations.

A second problem of convergence has been described in Sec. II: the difficulty to cross the barrier of the phase transition, particularly from liquid to solid. This transition involves the transformation of a collective order parameter. One way to enhance the efficiency of this transition would be to use precisely such collective and concerted motions in the MC exploration of the configurational space, 41,42 in view of simulating the real evolution of the system.

Finally, it has to be mentioned that we plan to explore the Fe system with the much more realistic MEAM potential, which takes into account all the environment of each atom. This potential being time consuming (three orders of magni-tude slower than the Rydberg one), an EAM type potential is under development. This approximate potential, only three times slower than the Rydberg one, will be used in the nested potential approach, 43,44 in which the fast EAM potential will alternate with the slower but much more sophisticated MEAM potential.

APPENDIX A: ANALYSIS OF THE LINEAR BIAS METHOD

Equation (11c) is assumed to be checked so that the histogram has a maximum E max . To optimize the determination of σ (E p ), the parameters of the bias must be chosen so that the Gaussian histogram H is approximately centered at the middle of the interval, at E α + α . The index α characterizing the interval for quantities such as the slope B are omitted below.

The condition that the Gaussian is well drawn on the interval means that

= K A 1/2 = K 2E B -σ 1/2 , ( A1
)
with the parameter K ≥ 1 (typically around 3-4). Furthermore, we constrain to be smaller than the value E of the interval center by a factor L , which is necessarily greater than 1, and must localize the histogram in a sufficiently welldefined zone of the σ curve, so that typically L 5 -10. These conditions give a lower bound for B

B -σ = 2K 2 L 2 E (A2a) = 2K 2 L 2 σ B T (A2b) ≥ 2L 2 σ B T , (A2c)
which is generally much smaller than the canonical value B T owing to the denominator σ . Then, this condition allows a much larger choice than the canonical bias. Now, the estimate of a maximum is a difficult task for a noisy curve. Then, it is made by cutting an arc at the top of the histogram, above exp(-K 2 M ), with K M being typically 1/2 -1/3. We develop ln[g/b L ] up to third order in energy around the maximum, taking into account a small asymmetry

H ( ) = H max × exp -A 2 + D 6 3 . (A3a) D = 2 B -σ E 2 + σ E . ( A3b
)
D is the third-order coefficient. This provides the expansion

H ( ) = H max × 1 -A 2 + D 6 3 + • • • . (A3c)
The contribution of the third-order D term will be shown as not negligible. The average of the energy over the delimited arc, E arc , is computed (this quantity being the least noised).

Finally, keeping only the lower terms in K 2 M , the maximum E max is expressed by

E max = E arc -D K 2 M 12A 2 . (A4)
Then, the corresponding value of the bias at this maximum is (A being expressed by Eq. (11b))

σ L b (E max ) = σ L b ( E arc ) + σ c1 + σ c2 , (A5a)
with the two contributions due to the third order,

σ c1 = - 2K 2 M 3 B B -σ , ( A5b
)

σ c2 = - K 2 M 3 Eσ B (B -σ ) 2 . ( A5c
)
When B is sufficiently greater than σ , σ c1 depends only on K M . We get -1/6 for K M = 1/2. The correction σ c2 is more involved and proportional to σ . The asymmetry D can be measured directly on the histogram arc by its third-order moment, but it is more accurate to compute an uncorrected σ curve, and from it to deduce the correction using Eqs. (A5b) and (A5c). We find that in some ranges, the correction σ c2 can reach 0.5-2, although it is generally much smaller.

Finally, an estimate of the error on E arc , due to the statistical fluctuations δH i of the histogram number H i , can be given. It will be made for simplicity linear least-square fits of the two sides of the arc defined by two values K Ma and K Mb . The accuracy can be expressed by

δ( E arc ) = 2 N arc dH dE i δH i , ( A6
)
where N arc is the number of bins of the arc, and dH /dE is the slope of the histogram at the limit of the arc. The evaluation of δH i can be made by remarking that the MC point wanders in a diffusion process from one side of the histogram to the other. Then, the counts accumulate around a bin i, leaves this region and returns to it after an average number of counts N RT . We assume that, for a total number of MC steps H T , δH i is proportional to the square root of the number of round-trips (ν RT = H T /N RT , N RT being the number of steps needed to go from one side to the other), and to the number of counts in the bin i during a round-trip (N RT H i /H T ) (A10) This result shows that δσ is minimized for B = 4σ , which is inside the condition of Eq. (11c). Reducing δE is also a way to improve the accuracy of δσ . Equation (A10) applied to the case E = 100 eV, δE = 0.1 eV, D E = 5 × 10 -5 eV 2 , K M = 1/2, H T = 10 9 counts, and B = 10 eV -1 gives δσ = 0.8, which is the correct magnitude of the observed fluctuations.

δH i = N RT H T 1/2 H i , ( A7

APPENDIX B: CORRECTION σ cd DUE TO DISCRETIZATION

For a bin of extension δE, centered in E, the histogram H is populated according to the ratio

H = δG δG b , ( B1
)
where δG = G(E + δE/2) -G(E -δE/2) is the number of states in the bin for the exact density of states (respectively δG b for the bias number of states). The correcting factor q g due to the non-zero extension δE can be defined as

δG = δE • g • q g , ( B2)
(respectively b and q b for the bias). Obviously, if δE → 0, q g → 1. With the previous two relations, the correction of σ due to the discretization can be derived

σ cd = -E d dE ln q g q b . ( B3
)
δG (respectively δG b ) can be expressed as δG = g(E)

+ δE 2 -δE 2 du exp E+u E dv σ g (v) v . ( B4
)
Approximating the exponential in terms of u by a constant integrant in terms of v, we get the result for q g in terms of the parameter x g q g = sh(x g ) x g , (B5a)

x g = σ g (E) δE 2E , ( B5b)
(respectively q b , x b , b(E)). This expression neglects higher order terms in 1/σ g . To avoid a too important correction, x 1 must be chosen, and in this case

q = sh(x) x = 1 + x 2 6 + • • • .
In this condition, using Eq. ( 13), the correction is expressed as

σ cd 2σ H (E) 3b(E) x b 2 . (B6)

FIG. 1 .

 1 FIG. 1. Histogram obtained with a σ -linear bias (σ L b (330) = 306, B = 40 eV -1 ) giving σ (349.93 eV) = 1103.20 for a 432 atom system.

FIG. 2 .

 2 FIG. 2. (a) Histograms obtained with the bias b (1) (respectively b (2) ) after the first (respectively second) iteration, for the range [230-1530 eV] above the phase transition. (b) Histograms obtained in the phase transition zone [150-220 eV] after 7.5 × 10 9 MC steps (1), 11.7 × 10 9 MC steps (2), 50 × 10 9 MC steps (3). All histograms are obtained for a 432 atom system.

FIG. 3 .

 3 FIG.3. Functions σ (E p ) (red curve) and β(E p ) (green curve, eV -1 ) for a 432 atom system, with two different energy scales: (a) from 0 to 400 eV and (b) from 400 to 6000 eV.

FIG. 5 .

 5 FIG. 5. Normalized curves: σ × 1000/N at as a function of E p × 1000/N at (eV) for 432, 1224, and 4044 atoms.

  5 and clarifies the meaning of the mapping function σ (E p ). When the equation β = β p (E p ), or equivalently σ (E p ) = β • E p , has a single intersection at the energy E β p , the integrant of Eq. (

FIG. 6 .

 6 FIG.6. Caloric curve U(T) (red line, eV), and heat capacity C V (T) (green line, in k B units) for a 432 atom system.

FIG. 7 .FIG. 8 .

 78 FIG. 7. Phase transition temperature T ph (K), in terms of N -1/3 at .

FIG. 9 .

 9 FIG. 9. Volumic density n(r) of a 432 atom cluster at different energies: 80, 135, 235, 650, 800, 950, 1100, and 1250 eV.

FIG. 10 .

 10 FIG.10. Pressure curve P(T) of a 432 atom cluster in semi-logarithmic scale. The left part depicts the zone of the triple point from 1000 to 3500 K. The right part is expanded up from 3500 to 120 000 K and shows the critical point.

FIG. 11 .

 11 FIG. 11. Histograms obtained in the phase transition zone [150-210 eV] for a 432 atom system. "This work" refers to the curve already presented in Fig. 2(b) with 50 × 10 9 MC steps. "WL + " refers to the result of a Metropolis calculation (100 × 10 9 MC steps) following a Wang-Landau process (27 × 10 9 MC steps).

FIG. 12 .

 12 FIG. 12. Evolution of the number of MC steps with ln f for a 432 atom system. Three energy ranges are investigated: [81, 171 eV] (solid), [150, 210 eV] (transition zone), and [231, 1530 eV] (liquid).

) 3 M 1 / 2 [

 312 Provided the system were out of the phase transition zone, it has been checked on various intervals E that N RT obeys approximately a diffusion process with a diffusion constant D E on the energy scale( E ) 2 = D E N RT . (A8)Choosing ( E ) 2 = 4/A, assuming in Eq. (A6) that the δH i are randomly distributed, and writing in Eq. (A7) that H i = H T • K M /N arc , the slope dH /dE at the limit of the arc is expressed asdH dE = 2 √ π δEH T K M A,δE being the width of the energy bin. Thus,δ( E arc ) arc A) 3/2 . (A9)Finally, N arc = 2K M /δE(A) 1/2 is eliminated and multiplied by B to obtainδσ ( E arc ) = πδE 2D E H T K 2E] 3/4 B [Bσ ] 3/4 .