使用周期性边界条件计算 Contact/Coordination 个数
Calculation of Contact/Coordination number with Periodic Boundary Conditions
我有 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 输出的不同，后者是用于 运行 模拟的软件。
感谢任何帮助。这里是否有一些意外行为，或者我错过了一些简单的事情？
对于这个例子，我发现：
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 接触。但我希望有 4 个连接
目前的理论是，当盒子足够薄以至于一个原子（具有周期性边界）可能与同一个原子接触不止一次时，就会发生此错误。现在正在寻找一种优雅的编码方式。我觉得这将只用 8 个原子解决这个问题。但是全模拟使用了2000个粒子和一个更大的盒子，还是有出入。
编辑：如下面的 P运行e 所示，情况确实如此，原子 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 添加校正 "skin" 因子并重新 运行 完整模拟 我发现我的值仍然经常低于真实结果，但在某些情况下它们会过冲.还是暗示有根本的区别？
是的，双击问题解释了 8 原子模拟的问题。双击对是 1-8 和 3-6，它们都通过 x 维度环绕。请注意其中一个维度的差异接近 Lx/2，而其他两个距离相对较小的实例。这意味着原子完全围绕你的 space-wrap 并接触到另一边。您当前的代码仅计算其中一次触摸。
要解决此问题，请查看 Lx 半径 - 半径范围内距离值的接触对（对 Ly 和 Lz 重复）。你找到的任何东西都需要 "inverted": distance = Lx-distance。然后你寻找另一组涉及仅你改变的触摸对。
请注意，只有当框尺寸小于半径的两倍时才会发生这种情况。你能为这种情况提供错误输出吗？由于您遇到了更大盒子的问题，因此上述问题不会解决更大套装的问题。也许有 8 个原子，但盒子大小为 2.1？
我有 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 输出的不同，后者是用于 运行 模拟的软件。
感谢任何帮助。这里是否有一些意外行为，或者我错过了一些简单的事情？
对于这个例子，我发现：
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 接触。但我希望有 4 个连接
目前的理论是，当盒子足够薄以至于一个原子（具有周期性边界）可能与同一个原子接触不止一次时，就会发生此错误。现在正在寻找一种优雅的编码方式。我觉得这将只用 8 个原子解决这个问题。但是全模拟使用了2000个粒子和一个更大的盒子，还是有出入。
编辑：如下面的 P运行e 所示，情况确实如此，原子 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 添加校正 "skin" 因子并重新 运行 完整模拟 我发现我的值仍然经常低于真实结果，但在某些情况下它们会过冲.还是暗示有根本的区别？
是的，双击问题解释了 8 原子模拟的问题。双击对是 1-8 和 3-6，它们都通过 x 维度环绕。请注意其中一个维度的差异接近 Lx/2，而其他两个距离相对较小的实例。这意味着原子完全围绕你的 space-wrap 并接触到另一边。您当前的代码仅计算其中一次触摸。
要解决此问题，请查看 Lx 半径 - 半径范围内距离值的接触对（对 Ly 和 Lz 重复）。你找到的任何东西都需要 "inverted": distance = Lx-distance。然后你寻找另一组涉及仅你改变的触摸对。
请注意，只有当框尺寸小于半径的两倍时才会发生这种情况。你能为这种情况提供错误输出吗？由于您遇到了更大盒子的问题，因此上述问题不会解决更大套装的问题。也许有 8 个原子，但盒子大小为 2.1？