把硬球装在盒子里

Question

我试着把硬球装在一个单位立方体的盒子里，这样这些球就不能相互重叠了。这是用Python实现的。 给我一些填充分数f，系统中的球数是N。 所以，我说每个球体的直径是 d = (p*6/(math.pi*N)**)1/3)。 我的盒子有周期性的边界条件——这意味着我的盒子在所有方向上都有一个重复出现的图像。如果有一个粒子位于长方体的边缘，并且有一部分粒子超出了墙，那么它将在另一侧突出

我的尝试：

1. 创建一个numpy N-by-3数组box，它保存每个粒子的位置向量[x,y,z]
2. 第一个粒子本来就很好
3. 阵列中的下一个粒子将与之前的所有粒子一起检查。如果它们之间的距离大于d，则继续移动到下一个粒子。如果它们重叠，随机更改所讨论粒子的位置向量。如果新位置没有与先前的原子重叠，则接受它
4. 对下一个粒子重复步骤2-3

我试图用这些硬球填充我的盒子，方式如下：

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

这段代码的问题是，如果p=0.45，则需要非常非常长的时间才能收敛。有没有更好的方法更有效地解决这个问题

Answer 1

我想你要找的是六边形密排晶格（HCP或有时称为面心立方晶格，FCC）或立方密排晶格（CCP）。例如，见Close-packing of equal spheres上的维基百科

由于你的空间有周期性条件，我相信你选择哪一个（hcp或ccp）并不重要，它们都达到了相同的密度约74.04%，这是高斯通过晶格堆积证明的最高密度

---

更新：

关于如何有效生成这样一个晶格的后续问题，让我们以HCP lattice为例。首先，让我们创建一组(i, j, k)索引[(0,0,0), (1,0,0), (2,0,0), ..., (0,1,0), ...]。然后，从这些索引中获取xyz坐标并返回数据帧：

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

我们可以使用plotly进行交互式探索，将结果绘制为scatter3d：

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()

注意：plotly的scatter3d不是一个很好的球体渲染：标记大小是恒定的（因此，当您放大和缩小时，“球体”将显示为改变相对大小），并且没有阴影、有限的z顺序忠实度等，但与绘图交互很方便

调整大小并剪辑到单元框中：

这里是一个严格的剪裁（每个球体都需要完全位于单位框内）。“周期性边界条件”是您需要单独解决的问题（请参阅下文了解更多想法）

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

这样，您会发现0.06的半径为608个完全封闭的球体：

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

你下一步要去哪里：

你可能会更深入地研究所谓的“周期性边界条件”的影响，也许会玩一些旋转（和小的平移）

为此，您可以尝试生成一个HCP晶格，该晶格足够大，任何旋转都将完全包围您的单位立方体。例如：

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

然后根据您的研究需要旋转（任意数量）和平移（任意方向最多1个半径），然后进行剪辑。“周期性边界条件”，正如您所说的，随着剪裁变得更加复杂，会带来一些额外的挑战。首先，剪辑中心位于长方体外部的任何球体。然后选择离边界足够近的球体，甚至沿立方体的墙将感兴趣的区域划分为重叠区域，然后检查每个区域中的球体之间的碰撞（根据周期性边界条件）