使用周期性边界条件计算接触/协调数

Question

我有 N 个原子的数据，包括每个原子的半径、每个原子的位置和盒子的尺寸。我想计算每个原子接触的原子数并存储结果。这相当于计算指定截止距离内的原子数。

盒子有周期性边界条件，因此在重新映射原子时会出现困难，以便计算出每个原子之间的距离尽可能小。

这是我目前所拥有的：

import numpy as np
Lx = 1.8609087768899999
Ly = 1.8609087768899999
Lz = 1.8609087768899999
natoms = 8
#Extract all atom positions
atomID = [1, 2, 3, 4, 5, 6, 7, 8]
x = [0.305153, 0.324049, 0.14708, 0.871106, 1.35602, 1.09239, 1.84683, 1.23874]
y = [1.59871, 1.40195, 0.637672, 0.859931, 1.68457, 0.651108, 0.694614, 1.57241]
z = [1.11666, 0.0698172, 1.60292, 0.787037, 0.381979, -0.00528695, 0.392633, 1.37821]
radius = np.array([0.5 for i in range(natoms)])

#Create matrix of particle distance
dx = np.subtract.outer(x,x)
dy = np.subtract.outer(y,y)
dz = np.subtract.outer(z,z)
radius_sum = np.add.outer(radius,radius)

#Account for periodic boundary conditions.
dx = dx - Lx*np.around(dx/Lx)
dy = dy - Ly*np.around(dy/Ly)
dz = dz - Lz*np.around(dz/Lz)
distance = np.array( np.sqrt( np.square(dx) + np.square(dy) + np.square(dz) ) )

touching = np.where( distance < radius_sum, 1.0, 0.0) #elementwise conditional operator (condition ? 1 : 0 in this case)
coordination_number = np.sum(touching, axis=0) - 1    # -1 to account for all diagonal terms being 1 in touching.

x、y、z 和 radius 是长度为 N（原子数）的数组。

此处计算的某些配位数与 LAMMPS 输出的配位数不同，LAMMPS 是用于运行仿真的软件。

感谢任何帮助。这里有一些意想不到的行为，还是我遗漏了一些简单的东西？

对于这个例子，我发现：

calculated    expected
3.0           4.0
4.0           4.0
3.0           4.0
4.0           4.0
4.0           4.0
4.0           5.0
5.0           5.0
3.0           4.0

我进一步计算了我的程序认为接触的原子，发现它们是：

1: [2, 4, 8]
2: [1, 3, 5, 7] 
3: [2, 6, 7]
4: [1, 6, 7, 8] 
5: [2, 6, 7, 8] 
6: [3, 4, 5, 7]
7: [2, 3, 4, 5, 6]
8: [1, 4, 5]

这是使用代码计算的：

contact_ID = [[] for i in range(natoms)]
for i in range(natoms):
    for j in range(natoms):
        if i==j: continue
        if (abs(touching[i,j] - radius_sum[i,j])) < .00001:
            contact_ID[atomID[i]-1].append(atomID[j])
        else: continue

我无法在模拟软件中做同样的事情。我可以看到 atomID 3 我的程序说 atom 3 与 atom 2、6 和 7 接触。但我预计会有 4 个连接

当前的理论是，当盒子足够薄以至于一个原子（具有周期性边界）可能与同一个原子多次接触时，就会发生此错误。现在正在寻找一种优雅的编码方式。我觉得这将解决这个问题，只有 8 个原子。但是完整的模拟使用了 2000 个粒子和一个更大的盒子，仍然存在差异。

编辑：如下面的 Prune 所示，情况确实如此，原子 3-6 和 1-8 实际上是双重接触的。伟大的！但是...这种情况已经非常不物理了，并且在我的模拟中永远不会发生，错误仍然会发生。

对于不应出现此问题的较大盒子尺寸，我的问题仍然存在。我可以重新创建问题的最简单的情况是 L=2.1 和 natoms=21。在这里我找到了配位数：

ID   calculated    expected
1   12.0    12.0
8   11.0    11.0
11  13.0    13.0
14  10.0    10.0
15  11.0    11.0
2   11.0    12.0
4   10.0    10.0
5   12.0    12.0
6   11.0    11.0
7   12.0    12.0
10  10.0    10.0
3   12.0    12.0
12  11.0    11.0
13  11.0    11.0
20  12.0    12.0
21  10.0    10.0
9   9.0 9.0
17  11.0    11.0
16  10.0    10.0
18  10.0    10.0
19  11.0    12.0

在这里，我们看到 atomID 为 2 和 19 的原子都有 11 个触点，而我应该计算它们有 12 个。两者都不足（并且是列表中唯一不正确的值）意味着原子 2 和 19 是接触。手动完成考虑边界条件的这些原子之间的最短距离：

import numpy as np
L=2.1
atomID = [2, 19]
x = [2.00551, 1.64908]
y = [0.146456, 1.90822]
z = [1.28094, 0.0518928]

#Manually checking these values we find:
dx = x[1] - x[0] #(= -0.35643)
dy = y[1] - y[0] #(= 1.761764)
dz = z[1] - z[0] #(= -1.2290472)

#Account for period BC's:
dx = dx          #(= -0.35643)
dy = dy - L      #(= -0.338236)
dz = dz + L      #(= 0.8709528)

distance = np.sqrt(dx**2 + dy**2 + dz**2)
print distance  #(1.000002358)

这意味着我的代码实际上没有错误......太棒了，模拟软件是否有理由将这些粒子视为接触，而实际上它们彼此相距 0.000002358？如何在不添加校正因子的情况下校正我的值？

为 radius_sum 添加校正“皮肤”因子并重新运行完整的模拟我发现我的值仍然经常低于真实结果，但在某些情况下它们会超过。仍然暗示存在根本区别？

Answer 1

是的，双击问题解释了 8 原子模拟的问题。双点触控对是 1-8 和 3-6，它们都包含在 x 维度上。请注意其中一个维度的差异接近 Lx/2，而其他两个距离相对较小的情况。这意味着原子完全围绕你的空间包裹并接触另一侧。您当前的代码只计算一次触摸。

要解决此问题，请查看 Lx 半径 - 半径范围内的距离值的接触对（对 Ly 和 Lz 重复）。您发现的任何内容都需要“反转”：距离 = Lx-距离。然后你寻找另一组涉及仅你改变的那些接触对的接触。

请注意，只有当盒子尺寸小于半径的两倍时，才会发生这种情况。你能为这种情况提供错误输出吗？由于您遇到了较大盒子的问题，因此上述问题无法解决较大盒子的问题。也许有 8 个原子，但盒子大小为 2.1？