Packing hard spheres in a box

user_14024882 0
Viewed 669

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.

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?

algorithm python
1 Answers

I think what you are looking for is either the hexagonal closed-packed (HCP or sometime called face-centered cubic, FCC) lattice or the cubic closed-packed one (CCP). See e.g. Wikipedia on Close-packing of equal spheres.

Since your space has periodic conditions, I believe it doesn't matter which one you chose (hcp or ccp), and they both achieve the same density of ~74.04%, which was proved by Gauss to be the highest density by lattice packing.

Update:

For the follow-up question on how to generate efficiently one such lattice, let's take as an example the HCP lattice. First, let's create a bunch of (i, j, k) indices [(0,0,0), (1,0,0), (2,0,0), ..., (0,1,0), ...]. Then, get xyz coordinates from those indices and return a DataFrame with them:

def hcp(n):
    dim = 3
    k, j, i = [v.flatten()
               for v in np.meshgrid(*([range(n)] * dim), indexing='ij')]
    df = pd.DataFrame({
        'x': 2 * i + (j + k) % 2,
        'y': np.sqrt(3) * (j + 1/3 * (k % 2)),
        'z': 2 * np.sqrt(6) / 3 * k,
    })
    return df

We can plot the result as scatter3d using plotly for interactive exploration:

import plotly.graph_objects as go

df = hcp(12)

fig = go.Figure(data=go.Scatter3d(
    x=df.x, y=df.y, z=df.z, mode='markers',
    marker=dict(size=df.x*0 + 30, symbol="circle", color=-df.z, opacity=1),
))
fig.show()

enter image description here

Note: plotly's scatter3d is not a very good rendering of spheres: the marker sizes are constant (so when you zoom in and out, the "spheres" will appear to change relative size), and there is no shading, limited z-ordering faithfulness, etc., but it's convenient to interact with the plot.

Resize and clip to the unit box:

Here, a strict clipping (each sphere needs to be completely inside the unit box). Your "periodic boundary condition" is something you will need to address separately (see further below for ideas).

def hcp_unitbox(r):
    n = int(np.ceil(1 / (np.sqrt(3) * r)))
    df = hcp(n) * r
    df += r
    df = df[(df <= 1 - r).all(axis=1)]
    
    return df

With this, you find that a radius of 0.06 gives you 608 fully enclosed spheres:

hcp_unitbox(.06).shape  # (608, 3)

Where you would go next:

You may dig deeper into the effect of your so-called "periodic boundary conditions", and perhaps play with some rotations (and small translations).

To do so, you may try to generate an HCP-lattice that is large enough that any rotation will still fully enclose your unit cube. For example:

r = 0.2  # example
n = int(np.ceil(2 / r))
df = hcp(n) * r - 1

Then rotate it (by any amount) and translate it (by up to 1 radius in any direction) as you wish for your research, and clip. The "periodic boundary conditions", as you call them, present a bit of extra challenge, as the clipping becomes trickier. First, clip any sphere whose center is outside your box. Then select spheres close enough to the boundaries, or even partition the regions of interest into overlapping regions along the walls of your cube, then check for collisions among the spheres (as per your periodic boundary conditions) that fall in each such region.

user_7581740
Related