# How to calculate displacement in each coordinate and wrap w.r.t. lattice parameter?

The problem is my box is not orthorhombic (base angles different than 90 degrees). So, I am getting wrong results with the formula,

def distance(a, b):
dx = abs(a[0] - b[0])
x = min(dx, abs(A - dx))

dy = abs(a[1] - b[1])
y = min(dy, abs(B - dy))

dz = abs(a[2] - b[2])
z = min(dz, abs(C - dz))

return sp.sqrt(x**2 + y**2 + z**2)


Any suggestion on how to get the proper distance between any atom and its neighbours, including the periodic images of the neighbours. Basically, it makes the results of RDF and CN unphysical.

• Have you tried using the Ovito tool? There are built in functions to evaluate the RDF and CN from the menu. Alternatively, it might be possible to use a tool such as Vesta and convert the unit cell to a rectangular one.
– PBH
Dec 1 '21 at 1:33
• @PBH why not write an answer? Comments are temporary and can be removed permanently without any notification. Dec 1 '21 at 16:48

Sticking with ase you can use the built-in get_distance function. It supports non-orthogonal unit cells.

An example to get the distance between atom number a and atom number b is:

import ase.io
distance_a_to_b = atom.get_distance(a, b)


Additionally if you want to compute the Radial Distribution Function (RDF) and Coordination Number (CN), asap3 is a python package building on top of ase that natively does it.

• I used the asap3 for the RDF calculation but the results from it doesn't give correct CN. Also, surprisingly, the RDF for the same file is not the same for OVITO and ASAP3, with same bin size, cutoff etc. Dec 2 '21 at 3:49
• Have you checked that the cell was properly read by ase (atom.get_cell()) and periodic boundary conditions set up (atom.get_pbc()) ?
– Hebo
Dec 2 '21 at 14:54
• Hi, @hebo, upon checking I found that the cell is read properly bu ASE along with the pbc. Even more baffling. Dec 14 '21 at 7:14

Have you tried using the Ovito where there are built in functions to evaluate both RDF and CN of a system?

Additionally, it might be possible to use Vesta to convert the unit cell to a rectangular one (see this video). That way you can use your code to do the analysis.