Hard spheres coupling.py

From Werner KRAUTH

import random

N = 4
L1 = [[0.2, 0.2], [0.7, 0.2], [0.2, 0.7], [0.7, 0.7]] # replace by 1000 valid
L2 = [[0.3, 0.3], [0.8, 0.3], [0.3, 0.8], [0.8, 0.8]] # random configurations
                                                      # of N particles
sigma = 0.1
sigma_sq = sigma ** 2
delta = 0.1
iter = 0
while True:
    iter += 1
    k = random.randint(0, N - 1)
    b = [random.uniform(delta, 1.0 - delta), random.uniform(delta, 1.0 - delta)]
    a = L1[k]
    min_dist = min((b[0] - c[0]) ** 2 + (b[1] - c[1]) ** 2 for c in L1 if c != a)
    box_cond = min(b[0], b[1]) < sigma or max(b[0], b[1]) > 1.0 - sigma
    if not (box_cond or min_dist < 4.0 * sigma ** 2):
        a[:] = b
    a = L2[k]
    min_dist = min((b[0] - c[0]) ** 2 + (b[1] - c[1]) ** 2 for c in L2 if c != a)
    box_cond = min(b[0], b[1]) < sigma or max(b[0], b[1]) > 1.0 - sigma
    if not (box_cond or min_dist < 4.0 * sigma ** 2):
        a[:] = b
    print(L1)
    print(L2)
    print()
    if L1 == L2:
        print(iter)
        break