# Krauth 2002

W. Krauth Disks on a Sphere and two-dimensional Glasses arXiv:cond-mat/0209391 (2002)

# Paper

Abstract I describe the classic circle-packing problem on a sphere, and the analytic and numerical approaches that have been used to study it. I then present a very simple Markov-chain Monte Carlo algorithm, which succeeds in finding the best solutions known today. The behavior of the algorithm is put into the context of the statistical physics of glasses.

# Context

The material in this paper is taken up in Section 7.1 of my 2006 book "Statistical Mechanics: Algorithms and Computations".

# Python implementation

The following Python program is from my 2015 MOOC on Statistical Mechanics

```import random, math

def unit_sphere():
x = [random.gauss(0.0, 1.0) for i in range(3)]
norm =  math.sqrt(sum(xk ** 2 for xk in x))
return [xk / norm for xk in x]

def minimum_distance(positions, N):
dists = [math.sqrt(sum((positions[k][j] - positions[l][j]) ** 2 \
for j in range(3))) for l in range(N) for k in range(l)]
return min(dists)

def resize_disks(positions, r, N, gamma):
Upsilon = minimum_distance(positions, N) / 2.0
r = r + gamma * (Upsilon - r)
return r

N = 12
r = 0.3
nsteps = 500000
sigma  = 0.25
gamma  = 0.15
while True:
positions = [unit_sphere() for j in range(N)]
if minimum_distance(positions, N) > 2.0 * r: break
n_acc = 0
for step in xrange(nsteps):
k = random.randint(0, N - 1)
newpos = [positions[k][j] + random.gauss(0, sigma) for j in range(3)]
norm = math.sqrt(sum(xk ** 2 for xk in newpos))
newpos = [xk / norm for xk in newpos]
new_min_dist = min([math.sqrt(sum((positions[l][j] - newpos[j]) ** 2 \
for j in range(3))) for l in range(k) + range(k + 1, N)])
if new_min_dist > 2.0 * r:
positions = positions[:k] + [newpos] + positions[k + 1:]
n_acc += 1
if step % 100 == 0:
acc_rate = n_acc / float(100)
n_acc = 0
if acc_rate < 0.5:
sigma *= 0.5
r = resize_disks(positions, r, N, gamma)
R = 1.0 / (1.0 / r - 1.0)
density = 1.0 * N / 2.0 * (1.0 - math.sqrt(1.0 - r ** 2))
print 'final r', r
print 'final R', R
print 'final density', density
```