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

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

#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 )
  • $\begingroup$ I would try to solve the problem for mass and spring interactions. Have you tried that? $\endgroup$
    – nicoguaro
    Commented Jun 8, 2018 at 3:04

1 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.

hexagonal lattice overlaid with rectangular and rhombic boxes

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.