How are you calculating neighbours in matlab ?. What algorithm are you using ?. The best way to check this is to download voro ++ source code compile and get the executable. The command line option is pretty easy to use.

Thank you Arun. I tried with voro++. It came out that for 2D when it calculates neighbouring lists of a particle it considers same particles up and down.
For example, in the output file
9939 5479 6085 5781 6085 5481 9585 10034
…
….
…
In the first line the Id of the particles is 6085 and in the neighbor list it appears twice where it should not be in the neighbour list at all. Consequently, it shows that it has 8 neighbours. I need the voronoi area. I am suspicious whether it will give the perfect area.
In matlab there is one function called voronoi which calculates number of neighbours.

Yes, I am using dimension 2 in LAMMPS. But when the simulation is run it ask for length of the box along z. So, I put a box length. But During run of the simulation there is no change of the z coordinates.

units lj

atom_style atomic

dimension 2

boundary p p p

read_data system1.data

and the starting of the system1.data file is

LAMMPS data file via write_data, version 15 May 2015, timestep = 10

This is just a feature of the tessellation method. Whenever a Voronoi cell extends to a boundary, a cell face is created at the boundary and this is treated as another neighbor. Hence in your 2d example, the upper and lower free surfaces are treated as two neighbors. I will update the documentation to make this clearer.

Dear Aidan,
Thank you very much. In addition to that I need voronoi cell area. According to the calculation, Voronoi cell area (in 2D) is changed with the length of the box along the z dimension which is not expected in case 2D. Is it like the value which is obtained as voronoi cell area, it is basically volume and we need to divide it by the box height to get the cell area?