Interatomic distance-periodic boundary conditions-non cubic unit cell
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
https://scicomp.stackexchange.com/questions/20165/periodic-boundary-conditions-for-triclinic-box
https://scicomp.stackexchange.com/questions/3107/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 -LYu2 elif currAtom[1]<LY: modPos = currAtom +LYu2 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. 509