具有周期边界条件的接触/协调数的计算
我有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和半径是长度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，我的程序说原子3与原子2，6和7有联系。
当前的理论是，当盒子足够薄，一个原子(具有周期边界)可能不止一次地与同一个原子接触时，就会发生这种错误。正在寻找一种优雅的方法来编写这个代码。我觉得这只需要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中添加一个修正的“皮肤”因子并重新运行完整的模拟，我发现我的值仍然常常低于真实的结果，但是在某些情况下它们过高了。还意味着有根本的区别吗？
回答 1
Stack Overflow用户
发布于 2017-02-22 18:27:12
是的，双触问题解释了8原子模拟的问题.双触对是1-8和3-6，它们都是通过x维数进行包装的.注意其中一个维度的差异接近Lx/2的实例，而其他两个距离相对较小。这意味着原子完全围绕着你的太空包裹，并在另一边接触。您当前的代码只计算其中一个触点。
要解决这一问题，请查看范围Lx-半径-半径范围内的距离值的接触对(对Ly和Lz重复)。任何你认为需要“倒置”的东西:距离=Lx-距离。然后，你寻找另一组接触，只涉及那些你改变了的接触对。
请注意，只有当方框尺寸小于半径的两倍时，才能发生这种情况。您能为这种情况提供错误输出吗？由于您得到的问题是一个更大的盒子，上面的问题不会解决问题与更大的一组。也许有8个原子，但盒子的大小是2.1？
https://stackoverflow.com/questions/42397870
复制相似问题
腾讯云开发者
Copyright © 2013 - 2026 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有
深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 
粤公网安备44030502008569号
腾讯云计算（北京）有限责任公司 京ICP证150476号 | 京ICP备11018762号