Is there a more efficient hard sphere packing algorithm?

Question

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:
            print("this is bad")
##########

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?

Answer 1

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.

What you can plug in your code right now to get an answer: https://pynndescent.readthedocs.io/en/latest/

Further reading: https://www.cs.princeton.edu/cass/papers/www11.pdf

p.s. This is probably the most universal and complex algorithm there is. From selecting next film to watch, to determening if fingerprint matches.

Answer 2

Using the Lubachevsky-Stillinger algorithm is the best choice. You get up to 64% space filling for a random monodisperse sphere packing. Another classical study is the algorithm by Jodrey and Tory (paper).

1. W. S. Jodrey and E. M. Tory Computer simulation of close random packing of equal spheres, Phys. Rev. A 32, 2347 (1985)

Answer 3

For those looking for random packing, I suggest you look at https://github.com/VasiliBaranov/packing-generation (I am not affiliated).

It supports the Lubachevsky–Stillinger, Jodrey–Tory, and force-biased generation algorithms; it can calculate the particle insertion probability, Steinhardt Q6 global and local order measures, coordination numbers for non-rattler particles, pair correlation function, structure factor, and reduced pressure after pressure equilibration. It doesn't require any preinstalled libraries and is multiplatform (Windows/nix).

In particular, there is the -lsgd flag:

-lsgd: Lubachevsky–Stillinger with gradual densification: the LS generation is run until the non-equilibrium reduced pressure is high enough (e.g., a conventional value of 1e12), then compression rate is decreased (devided by 2) and the LS generation is run again, until the pressure is high enough again; this procedure is repeated until the compression rate is low enough (1e-4). See Baranau and Tallarek (2014) Random-close packing limits for monodisperse and polydisperse hard spheres, doi:10.1039/C3SM52959B.

You can specify an arbitrary particle size distribution by specifying the diameter of each particle, and you can also specify the initial locations of the particles if you prefer that over random location assignment. The code takes some fiddling with the diameters for post-processing, and it seems best in terms of time/performance tradeoff if you run certain algorithms in succession (-fba, then -ls, then -lsgd).

I put together a Jupyter notebook that should run a random hard sphere packing simulation to completion.