12
$\begingroup$

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?

$\endgroup$
14
  • 1
    $\begingroup$ 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) $\endgroup$ Dec 8 '20 at 4:58
  • $\begingroup$ 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? $\endgroup$
    – megamence
    Dec 8 '20 at 5:06
  • 1
    $\begingroup$ 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? $\endgroup$ Dec 8 '20 at 5:44
  • $\begingroup$ @TristanMaxson yes, it is a cubic box, with X, Y, Z going from 0 to 1 $\endgroup$
    – megamence
    Dec 8 '20 at 6:56
  • $\begingroup$ Final request, can you supply something executable that I can start with? $\endgroup$ Dec 8 '20 at 14:38
4
$\begingroup$

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: https://hexnet.org/files/images/hexnet/content_close-packing.png 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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.