Event-driven Monte Carlo algorithm for general potentials

In the present paper, we generalize the event-chain algorithm to the case of general potentials, greatly extending the scope of the original method. As for eventdriven molecular dynamics (MD), we discretize the potential (here on an energy scale ∆ E ), but unlike the MD algorithm [3,4], simulations at arbitrary small ∆ E are feasible.

ℓ Figure 1: Event-chain algorithm for hard disks [1]. Left : Collective move, with a total displacement ℓ as a fixed parameter of the algorithm. The move is never rejected. Right : Return move, demonstrating microreversibility.

In the event-chain algorithm [1], a randomly chosen disk is moved in a straight line by a distance ℓ, if this generates no overlap for all ℓ ′ ≤ ℓ, or else until it hits another particle i ′ . The latter particle, i ′ , then moves for the remainder of the distance ℓ or, again, until it hits yet another particle, etc. (see Fig. 1). The event-chain move is micro-reversible, and it satisfies detailed balance if all the displacements are in the same direction. For hard disks, it is about two orders of magnitude faster than the local Monte Carlo algorithm [5], and it compares favorably with the MD method [6].

In the event-driven Monte Carlo algorithm, we consider a discretized potential V (r), with discontinuities such that V (r) and thus the total potential energy are multiples of a given energy step ∆ E [7]. From an initial configuration of energy E, a randomly chosen particle i is moved in a straight line until the total displacement, at energy E only, equals a fixed value ℓ or else until it “hits” another particle i ′ . This means that the energy of the system would exceed E if the particle i moved any farther. The latter particle, i ′ , then moves for the remainder of the distance ℓ or, again, until it hits yet another particle, etc. (see Fig. 2). In the absence of boundary effects, the move is guaranteed to terminate. It is again rejectionfree and micro-reversible. Detailed balance may again be broken and the algorithm can be run with moves, say, in the +x and +y direction only. Particle i moves from its initial configuration until its collision event with particle i ′ (potential of Eq. ( 1) with ∆E = 1/4 and E = 3/4). Only displacements at energy E are discounted from ℓ.

The event-driven Monte Carlo algorithm can be im-plemented by computing the intersection points of the potential discontinuities with the trajectory of the particle i. One recovers E(x) by adding up the signs in the sorted list of the intersections. This is illustrated in Fig. 2 for the potential defined in Eq. ( 1). In Fig. 3, we consider particles interacting with the potential

for ∆ E = 1/10. We compare the pair-correlation function of the event-driven Monte Carlo algorithm for a large system with N = 128 2 particles to the results obtained with a local Monte Carlo algorithm in the same system. With the local Monte Carlo algorithm, only configurations with potential energy E are used. There is no doubt about the correctness of the algorithm and the implementation. With this simple implementation, we reached ∼ 5 × 10 9 collisions per hour on a 3 GHz workstation. In the naive implementation suggested in Fig. 2, the number of steps between events scales ∝ 1/∆ E , and the algorithm becomes very slow in the limit ∆ E → 0. Unlike for event-driven molecular dynamics [3], this difficulty can be overcome, and the speed of the algorithm remains constant for ∆ E → 0. The algorithm is straightforward for a convex potential with finite support, as for example the continuum limit of Eq. ( 1), given by

for which E(x) is continuous and piecewise C ∞ while ∂E/∂x is piecewise monotonously decreasing (see Fig. 4).

The root E cont (x root ) = E of the continuous potential is uniquely determined via a decision tree in subsequent C ∞ intervals of E(x). For finite ∆ E , one must only identify the discontinuity steps of the potential in the interval x root ± const[∂E/∂x(x root )] -1 . We have implemented the algorithm for very large values of 1/∆ E (with ℓ ∝ ∆ E ), and achieved constant scaling of the algorithm. Although the event-driven Monte Carlo algorithm can be implemented for arbitrarily small ∆ E , and does not

Figure 4: Energy E(x) for the configuration of Fig. 2, for different discretizations, and in the continuum limit, and derivative ∂Econt/∂x which is piecewise monotonously decreasing.

slow down then, it is interesting to consider the ∆ E → 0 limit. A good model for this is to suppose that rather than being fixed, ∆ -1 E is sampled, between events, from a distribution with very large mean. The collision rules in Fig. 2 then become probabilistic. They are readily analyzed in simple cases. The event itself, which takes place at energy E, is initiated and ended by particles coming in from and moving out to lower energies. This event may involve a sizable number of particles and thus become quite complex. It is unclear whether truly general stochastic rules for the collisions and the use of ℓ can be obtained.

In conclusion, in this paper, we have generalized the successful event-chain algorithm from hard spheres to the case of general potentials and provided details on its implementation. This algorithm is an alternative to the local Monte Carlo algorithm [5] and, even with a simple implementation, should outperform it in most applications. Moreover, arbitrarily small discretization steps can be handled efficiently, which is not the case in the eventdriven MD [3,4,6].