Dress Krauth 1995
From Werner KRAUTH
C. Dress, W. Krauth Cluster Algorithm for hard spheres and related systems Journal of Physics A: Math. Gen. 28, L597 (1995)
Abstract: In this paper, we present a cluster algorithm for the simulation of hard spheres and related systems. In this algorithm, a copy of the configuration is rotated with respect to a randomly chosen pivot point. The two systems are then superposed, and clusters of overlapping spheres in the joint system are isolated. Each of these clusters can be flipped independently, a process which generates non-local moves in the original configuration. A generalization of this algorithm (which works perfectly well at small density) can be successfully made to work at densities around the solid-liquid transition point in the two-dimensional hard-sphere system.
Further insights
In the original version, we formulated the algorithm as an exchange between two copies of a configuration. The "pocket" version of the algorithm, that I also used, with R. Moessner, in a 2003 paper on dimer models, is much easier. This pocket version that is implemented in the below Python program.
Python version of the cluster algorithm
###========+=========+=========+=========+=========+=========+=========+= ## PROGRAM : pocket_disks.py ## PURPOSE : implement SMAC algorithm 2.18 (hard-sphere-cluster) ## for four disks in a square of sides 1 ## LANGUAGE: Python 2.5 ##========+=========+=========+=========+=========+=========+=========+ from random import uniform as ran, choice def box_it(x): x[0]=x[0]%1. if x[0] < 0 : x[0]=x[0]+1 x[1]=x[1]%1. if x[1] < 0 : x[1]=x[1]+1 return x def dist(x,y): d_x= abs(x[0]-y[0])%1 d_x = min(d_x,1-d_x) d_y= abs(x[1]-y[1])%1 d_y = min(d_y,1-d_y) return d_x**2 + d_y**2 def T(x,Pivot): x=(2*Pivot[0]-x[0],2*Pivot[1]-x[1]) return x # # Program starts here # Others=[(0.25,0.25),(0.25,0.75),(0.75,0.25),(0.75,0.75)] sigma_sq=0.15**2 for iter in range(10000): a = choice(Others) Others.remove(a) Pocket = [a] Pivot=(ran(0,1),ran(0,1)) while Pocket != []: a = choice(Pocket) Pocket.remove(a) a = T(a,Pivot) for b in Others[:]: # "Others[:]" is a copy of "Others" if dist(a,b) < 4*sigma_sq: Others.remove(b) Pocket.append(b) Others.append(a) print Others, ' ending config '