# Is there a more efficient hard sphere packing algorithm?

I am trying to pack hard-spheres in a unit cubical box, such that these spheres cannot overlap on each other. This is being done in Python. I am given some packing fraction f, and the number of spheres in the system is N. So, I say that the diameter of each sphere will be d = (p*6/(math.pi*N)**)1/3). My box has periodic boundary conditions - which means that there is a recurring image of my box in all direction. If there is a particle who is at the edge of the box and has a portion of it going beyond the wall, it will stick out at the other side.

This is my code, with me generating a configuration array called box.

import numpy as np
import pandas as pd
import math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

np.random.seed(1)
#initialize the system
L = 1 #length of box
N = 500 #number of particles in a box
phi = np.asarray([0.2,0.45]) #packing fractions
#calculating the particle diameters so that we get the right density and packing fractions
d = (phi[1]*6/(N*math.pi))**(1/3)
#d=0.12 (approx)
rho = N/(L**3) #number density

#intitalize box
box = np.random.uniform(0,1,(N, 3))

#simulation 1
diameter = d #tolerance
#check if box is valid

#this is a collection of loops to check if the particles we have placed are
#not overlapping
#we place i = 0 particle, no problem
#then we place the other particles, and we check with all the previously placed particles if
#the placement is okay or not

for i in range(1,N):
mybool=True
print("particles in box: " + str(i))
while (mybool): #the deal with this while loop is that if we place a bad particle, we need to change its position, and restart the process of checking
for j in range(0,i):
displacement=box[j,:]-box[i,:]
for k in range(3):
if abs(displacement[k])>L/2:
displacement[k] -= L*np.sign(displacement[k])
distance = np.linalg.norm(displacement,2) #check distance between ith particle and the trailing j particles
if distance<diameter:
box[i,:] = np.random.uniform(0,1,(1,3)) #change the position of the ith particle randomly, restart the process
break
if j==i-1 and distance>diameter:
mybool = False
break
#print(np.linalg.norm(box[0,:]-box[119,:]))
#test script to check if the above generated a sound configuration worked or not
for i in range(1,N):
for j in range(0,i):
displacement=box[j,:]-box[i,:]
for k in range(3):
if abs(displacement[k])>L/2:
displacement[k] -= L*np.sign(displacement[k])
distance = np.linalg.norm(displacement,2) #check distance between ith particle and the trailing j particles
if distance<diameter:
##########


My attempt:

1. Create a numpy N-by-3 array box which holds the position vector of each particle [x,y,z]
2. The first particle is fine as it is.
3. The next particle in the array is checked with all the previous particles. If the distance between them is more than d, move on to the next particle. If they overlap, randomly change the position vector of the particle in question. If the new position does not overlap with the previous atoms, accept it.
4. Repeat steps 2-3 for the next particle.

I am trying to populate my box with these hard spheres, in the following manner:

for i in range(1,N):
mybool=True
print("particles in box: " + str(i))
while (mybool): #the deal with this while loop is that if we place a bad particle, we need to change its position, and restart the process of checking
for j in range(0,i):
displacement=box[j,:]-box[i,:]
for k in range(3):
if abs(displacement[k])>L/2:
displacement[k] -= L*np.sign(displacement[k])
distance = np.linalg.norm(displacement,2) #check distance between ith particle and the trailing j particles
if distance<diameter:
box[i,:] = np.random.uniform(0,1,(1,3)) #change the position of the ith particle randomly, restart the process
break
if j==i-1 and distance>diameter:
mybool = False
break


The problem with this code is that if p=0.45, it is taking a really, really long time to converge. Is there a better method to solve this problem, more efficiently?

• Do you need this box to be packed randomly? You might fill it in some more optimal solution then perturb the entire box (think high temperature MD of elastic spheres) – Tristan Maxson Dec 8 '20 at 4:58
• no the box need not be packed randomly. So are you saying I should arrange my spheres in a lattice and perturb them? What is the best way to arrange spheres in an optimal way and perturbing them? – megamence Dec 8 '20 at 5:06
• My initial response is that at the very worse, you can FCC pack a cube into the corner then fill randomly from there. Not suggesting its a good idea, but it would be a much faster solution. Is your goal just to get a faster solution within python? This is technically somewhat of a stack overflow question if so, but I think we can handle it. Wait is this a cubic box? X,Y,Z are all orthogonal and the same length? – Tristan Maxson Dec 8 '20 at 5:44
• @TristanMaxson yes, it is a cubic box, with X, Y, Z going from 0 to 1 – megamence Dec 8 '20 at 6:56
• Final request, can you supply something executable that I can start with? – Tristan Maxson Dec 8 '20 at 14:38