Interatomic distance-periodic boundary conditions-non cubic unit cell

Question

I am trying to find interatomic distance considering periodic boundary conditions for hexagon cubic cells (graphite). I tried to follow the answers to these two questions here but am unable to get the right results

periodic boundary conditions for triclinic box

Minimum image convention for triclinic unit cell

My code, implementing the answer at first post is 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 )

The code for the second post is

#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 )

Answer 1

I think that you may be overcomplicating things, trying to use general triclinic periodic boundary conditions for a system with 2D hexagonal periodicity. It is quite possible to use a rhombic cell in the $x$ and $y$ directions (with regular boundaries in $z$ if needed, making it a rhombic prism in 3D) or indeed a rectangular cell (i.e. a cuboidal box in 3D). Both these are illustrated in the figure.

For the blue box illustrated here, if the side of the hexagon is one unit, then $L_x=12$ and $L_y=6\sqrt{3}\approx 10.392$, but there is plenty of flexibility depending on how many unit cells (in each direction) you wish to include. Then the Python code for periodic boundary corrections is the same as for any cuboidal box, for instance xij=xij-np.rint(xij/Lx)*Lx and similarly for yij and zij, where (xij,yij,zij) is an uncorrected pair vector.