问题描述
我有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接触.但是我希望有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
在这里，我们看到原子ID为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中添加一个校正皮肤"因子并重新运行完整的模拟，我发现我的值仍然经常低于真实结果，但是在某些情况下它们会超出实际值.仍然暗示存在根本差异吗?
是的，两次接触问题解释了8原子模拟的问题.双触对是1-8和3-6，它们都通过 x 维包裹.注意，其中一个维度的差异接近Lx/2，而其他两个距离相对较小的情况.这意味着原子完全到达您的空间包裹物周围，并在另一侧接触.您当前的代码只算一次接触.
要解决此问题，请查看触摸对以获取Lx半径-半径范围内的距离值(重复Ly和Lz).您发现的任何东西都需要反转":distance = Lx-distance.然后，您查找另一组涉及您仅更改过的触摸对的触摸.
请注意，只有当框尺寸小于半径的两倍时，才会发生这种情况.在这种情况下，您可以提供错误输出吗?由于您遇到的是较大包装盒的问题，因此上述问题不会解决较大包装盒的问题.也许有8个原子的东西，但是盒子的大小是2.1?
I have data for N atoms including the radius of each atom, the position of each atom and the box dimensions. I want to calculate the number of atoms each atom is in contact with and store the result. This is equivalent to calculating the number of atoms within a specified cutoff distance.
The box has periodic boundary conditions, and hence a difficulty arises in remapping atoms such that the distance calculated between each atom is the minimum possible.
Here is what I have so far:
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 and radius are arrays of length N (number of atoms).
Some of the coordination numbers calculated here differ to the ones LAMMPS outputs, which is the software used to run the simulation.
Any help is appreciated. Is there some unexpected behaviour here, or am I missing something simple?
For this example I find:
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
I have further calculated the atoms which my program thinks are in contact and found them to be:
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]
This was calculated using the code:
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
There is no way for me to do the same from within the simulation software. I can see for atomID 3 my program says atom 3 is in contact with atoms 2, 6 and 7. But I expect there to be 4 connections
Current theory is that this error is occurring when the box is thin enough that one atom (with periodic boundaries) could be in contact with the same atom more than once. Looking for an elegant way to code this right now. I feel this will fix this problem with only 8 atoms. However the full simulation uses 2000 particles and a much larger box, and there still exist discrepancies.
EDIT: As shown by Prune below this was indeed the case, atoms 3-6 and 1-8 were in fact in double contact. Great! However... this case is already highly unphysical and will never occur in my simulation, where the errors still occur.
My issue still remains for larger box sizes where this problem shouldn't occur. The most simple case I can recreate the problem for is with L=2.1 and natoms=21. Here I find coordination numbers of:
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
Here we see that atoms with atomID 2 and 19 both have 11 contacts, whereas I should have calculated them to have 12. Both falling one short (and being the only incorrect values in the list) imply that atoms 2 and 19 are in contact. Manually completing the shortest distance between these atoms accounting for boundary conditions:
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)
This implies there is in fact not an error in my code... Fantastic, is there a reason that the simulation software would count these particles as touching when in fact they are 0.000002358 away from each other? How am I to go about correcting my values without adding a correction factor?
Adding a correction "skin" factor to the radius_sum and rerunning the full simulation I find that my values still often undershoot the true result, however in some cases they overshoot. Still implies that there is a fundamental difference?
Yes, the double-touch issue explains the problem with the 8-atom simulation. Double-touch pairs are 1-8 and 3-6, both of those wrapping through the x dimension. Note the instances in which a difference in one of the dimensions is close to Lx/2, while the other two distances are relatively small. This means that the atoms reach fully around your space-wrap and touch on the other side. Your current code counts only one of the touches.
To fix this, look at touching pairs for distance values in the range Lx-radius - radius (repeating for Ly and Lz). Any you find need to be "inverted": distance = Lx-distance. Then you look for another set of touches involving only those touching pairs you changed.
Note that this can happen only when a box dimension is less than twice the radius. Can you provide error output for such a case? Since you're getting the problem for a larger box, the above issue won't solve the problem with the larger set. Perhaps something with 8 atoms, but a box size of 2.1?
这篇关于具有周期边界条件的接触/协调数的计算的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！