# 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*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) Dec 8, 2020 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? Dec 8, 2020 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? Dec 8, 2020 at 5:44
• @TristanMaxson yes, it is a cubic box, with X, Y, Z going from 0 to 1 Dec 8, 2020 at 6:56
• @megamence Why not answer the question with what you found then? You'll help others that might have the same question, in the years and decades and generations to come! Also our answer ratio needs improvement, which means we need to get more answers/question. Dec 17, 2020 at 1:47

If you want faster solution, you can use simplification. If there are more than some critical amount of spheres (about 100), then grid solution is likely almost as good as true solution. Make a triangle grid, and only check the grid size, adding or removing rows and columns. Add another grid on top, until you have enough spheres to satisfy your requirements. Most closely packed spheres pattern looks like this: And there is no need to calculate spheres in between the edges, because their perfect arrangement is known. You will have gaps on the borders of a box, even with 'poking from another side', because spheres size likely wont arrange to a perfect grid on the edge. This solution takes about the same time for any size of the box and number of spheres, because you are only checking the boundary.

If you want more precise solution, you need physical simulation. Put all your spheres randomly in a box, with overlapping. Then start the simulation with gravity and pushing overlapping spheres away. Calculate each step all the forced that act on a sphere, sum it, apply the sum as a force acting against sphere's inertia. You will need nearest neighbor search, that will be your true enemy, simulation itself is simple. https://en.m.wikipedia.org/wiki/Nearest_neighbor_search There are a few dozen solutions, but since you will need dynamic, lots of them wont work fast enough, mostly because tree rebuilding is often expensive. I would suggest graph based, local, without global index.