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?
