考虑到六边形立方晶胞（石墨）的周期性边界条件，我试图找到原子间距离。我试图在这里遵循这两个问题的答案，但无法获得正确的结果
我的代码，在第一篇文章中实现答案是here
positions = initialConfig.get_positions()
nAtoms = positions.shape[0]
xMin = zMin = 0.0
xMax = LX
zMax = LZ
yCoord = positions[:,0:2]
u2 = UCell[1][0:2,]
modYPos = np.zeros((0,2))
for i in range(nAtoms):
currAtom = yCoord[i,:]
if currAtom[1]>LY:
modPos = currAtom -LY*u2
elif currAtom[1]<LY:
modPos = currAtom +LY*u2
else:
modPos = currAtom
modYPos = np.vstack((modYPos, modPos))
xCoord = modYPos[:,0].reshape(-1,1)
xCoord = xCoord - xMax*np.rint(xCoord/xMax)
zCoord = positions[:,2].reshape(-1,1)
zCoord = zCoord - zMax*np.rint(zCoord/zMax)
app5rMat = np.zeros((nAtoms,coordinates))
app5rMat[:,0] = xCoord[:,0]
app5rMat[:,1] = yCoord[:,1]
app5rMat[:,2] = zCoord[:,0]
app5RIJ = app5rMat.reshape(nAtoms,1,coordinates)-app5rMat.reshape(1,nAtoms,coordinates)
app5Dist = np.linalg.norm(app5RIJ, axis =2 )
第二篇文章的代码是
#S1 create A matrix
UCell = initialConfig.get_cell()
cellParams = initialConfig.get_cell_lengths_and_angles()
a = cellParams[0]
b = cellParams[1]
c = cellParams[2]
alpha= math.radians(90)#math.pi/2#math.radians(cellParams[3])
beta= math.radians(90)#math.pi/3#math.radians(cellParams[4])
gamma= math.radians(60)#math.pi/2#math.radians(cellParams[5])
A = UCell.T
#S2 Invert A to bet B
B = np.linalg.inv(A)
np.matmul(B,A)
#S3 -- get fractional coordinates
positions = initialConfig.get_positions()
nAtoms = positions.shape[0]
fracCood = np.zeros((0,3))
for i in range(nAtoms):
currAtom = positions[i,:].reshape(-1,1)
fracCood = np.vstack((fracCood, np.dot(B,currAtom).T))
fracCood.shape
#s4 -- translate into reference cell
onesArray = np.ones((nAtoms,coordinates))
g = fracCood - np.rint(2*fracCood - onesArray)#np.floor(fracCood)
#s4-- translate into real space
realCood = np.zeros((0,3))
for i in range(nAtoms):
currAtom = g[i,:].reshape(-1,1)
realCood = np.vstack((realCood, np.dot(A,currAtom).T))
algoRij = g.reshape(nAtoms,1,coordinates)-g.reshape(1,nAtoms,coordinates)
secondPart = np.rint(2*algoRij.copy() - onesArray.reshape(nAtoms,1,coordinates))
algoRij = algoRij-secondPart
algoRijPer = algoRij #-np.rint(algoRij)
algoRijDis = np.linalg.norm(algoRijPer, axis =2 )
