From Werner KRAUTH

This page presents the program event_chain.py, a Markov-chain algorithm, naively implemented for four disks in a square box of sides 1 with periodic boundary conditions.

Reference

E. P. Bernard, W. Krauth, D. B. Wilson Event-chain algorithms for hard-sphere systems Physical Review E 80 056704 (2009)

Description

This is a Python implementation of the hard-disk event-disk algorithm, for four particles in a square box of length 1, with periodic boundary conditions. The algorithm belongs to a class of irreversible infinitesimal Markov chains using the lifting paradigm, as discussed in this article with Manon Michel and Sebastian Kapfer. In this article, we also discuss the lifting paradigm with infinitesimal moves also for interacting particles.

The event-chain algorithm was applied, very successfully, to the problem of hard-disk melting, in a 2011 paper with Etienne Bernard and to three-dimensional hard spheres, with Masaharu Isobe from the Nagoya Institute of Technology, in Japan.

Program

import random, math

def event(a, b, dirc, sigma):
    d_perp = abs(b[not dirc] - a[not dirc]) % 1.0
    d_perp = min(d_perp, 1.0 - d_perp)
    if d_perp > 2.0 * sigma:
        return float("inf")
    else:
        d_para = math.sqrt(4.0 * sigma ** 2 - d_perp ** 2)
        return (b[dirc] - a[dirc] - d_para + 1.0) % 1.0

L = [[0.25, 0.25], [0.25, 0.75], [0.75, 0.25], [0.75, 0.75]]
ltilde = 0.819284; sigma = 0.15
for iter in xrange(20000):
    dirc = random.randint(0, 1)
    print iter, dirc, L
    distance_to_go = ltilde
    next_a = random.choice(L)
    while distance_to_go > 0.0:
        a = next_a
        event_min = distance_to_go
        for b in [x for x  in L if x != a]:
            event_b = event(a, b, dirc, sigma)
            if event_b < event_min:
                next_a = b
                event_min = event_b
        a[dirc] = (a[dirc] + event_min) % 1.0
        distance_to_go -= event_min

Version

See history for version information.