Is the Neighbors List correct ?

[Optimox]

Referring to the discussion on #839 I am looking at voronoi bonds given by voro.nlist.

I'm interested in finding all the points which have a bond crossing on of the frontier of the box using:

unwrapped_bonds = centered_points[nlist.query_point_indices] \
                        + box.wrap(centered_points[nlist.point_indices] \
                        - centered_points[nlist.query_point_indices])
box.contains(unwrapped_bonds)

I then plot in red the bonds that cross the frontier and in blue the one that don't.

It seems that the method to differentiate between red and blue bonds work well. However I feel like some red bonds are missing:

Do you agree with me on some red bonds missing? Or am I not understanding something?

How exactly is the neighbors list computed? Is it supposed to be exhaustive?

Thanks in advance for your help

[tommy-waltmann]

Hi @Optimox ,

Can you give a script which details the box and points you are working with so we can look at this more in-depth?

Based on the discussion in #839 , you are working with a 3D box, but I think displaying it on a 2D rectangle can create distortions which may be misleading as to which cell is a neighbor of which.

[Optimox]

Hello @tommy-waltmann,

Thanks for your help.

Indeed my original problem is 3D, however in order check my pipeline I started with 2D. So I work with a tiny slice of my 3D data and put everything in 2D before doing any computation. So there are no 3D -> 2D artefacts in the drawing above.

You'll find bellow a complete code snippet which allow you to reproduce this exact drawing. It's a bit long sorry.

I'm really curious to understand what's happening, I hope you'll be able to help me.

import numpy as np
import freud
from matplotlib import pyplot as plt
import matplotlib


def set_size(w,h, ax=None):
    """ plot utility to have bigger image"""
    if not ax: ax=plt.gca()
    l = ax.figure.subplotpars.left
    r = ax.figure.subplotpars.right
    t = ax.figure.subplotpars.top
    b = ax.figure.subplotpars.bottom
    figw = float(w)/(r-l)
    figh = float(h)/(t-b)
    ax.figure.set_size_inches(figw, figh)
    
def check_inside_box(points, freud_box):
    """
        This is a simple numpy check that points belong to original cubic box.        
    """
    
    Lx, Ly, Lz = freud_box.Lx, freud_box.Ly, freud_box.Lz
    if Lz == 0:
        # for 2D don't take Lz into account
        Lz = 1000
    is_x = (points[:, 0] <= Lx/2) & (points[:, 0] >= -Lx/2)
    is_y = (points[:, 1] <= Ly/2) & (points[:, 1] >= -Ly/2)
    is_z = (points[:, 2] <= Lz/2) & (points[:, 2] >= -Lz/2)
    return (is_x & is_y & is_z)

centered_points = np.array([[-250,   22,    0],
       [-248,   14,    0],
       [-246,   24,    0],
       [-240,   28,    0],
       [-234,    7,    0],
       [-232,   17,    0],
       [-230,   11,    0],
       [-226,    5,    0],
       [-224,   20,    0],
       [-222,   -2,    0],
       [-220,   -1,    0],
       [-214,    5,    0],
       [-214,   17,    0],
       [-212,   11,    0],
       [-208,   10,    0],
       [-202,    4,    0],
       [-202,   28,    0],
       [-198,    5,    0],
       [-196,    9,    0],
       [-196,   16,    0],
       [-190,   21,    0],
       [-184,  -11,    0],
       [-182,   11,    0],
       [-178,  -30,    0],
       [-178,    8,    0],
       [-178,   22,    0],
       [-176,   -1,    0],
       [-174,   17,    0],
       [-168,  -12,    0],
       [-166,   11,    0],
       [-164,    2,    0],
       [-158,   25,    0],
       [-154,    6,    0],
       [-152,    1,    0],
       [-146,    9,    0],
       [-140,    4,    0],
       [-140,   14,    0],
       [-132,   -1,    0],
       [-130,   22,    0],
       [-124,    3,    0],
       [-120,   -5,    0],
       [-120,   10,    0],
       [-108,    9,    0],
       [-104,   10,    0],
       [ -96,   18,    0],
       [ -88,   -6,    0],
       [ -86,   -6,    0],
       [ -82,    7,    0],
       [ -78,  -13,    0],
       [ -78,    1,    0],
       [ -76,   13,    0],
       [ -74,   -4,    0],
       [ -66,   -2,    0],
       [ -62,   23,    0],
       [ -56,   -7,    0],
       [ -56,   13,    0],
       [ -54,    7,    0],
       [ -50,   16,    0],
       [ -50,   28,    0],
       [ -48,    0,    0],
       [ -46,  -19,    0],
       [ -46,    6,    0],
       [ -46,   21,    0],
       [ -42,    1,    0],
       [ -38,    7,    0],
       [ -36,   -3,    0],
       [ -36,   12,    0],
       [ -36,   26,    0],
       [ -32,    0,    0],
       [ -30,   14,    0],
       [ -30,   20,    0],
       [ -26,    5,    0],
       [ -24,    9,    0],
       [ -22,  -12,    0],
       [ -22,    2,    0],
       [ -20,   21,    0],
       [ -16,    8,    0],
       [ -14,   -6,    0],
       [ -14,    4,    0],
       [ -14,   16,    0],
       [ -10,    5,    0],
       [ -10,   28,    0],
       [  -8,    5,    0],
       [  -4,   -9,    0],
       [  -2,    5,    0],
       [  -2,   15,    0],
       [   0,   -4,    0],
       [   0,   22,    0],
       [   4,   17,    0],
       [   6,    2,    0],
       [  10,    1,    0],
       [  12,   26,    0],
       [  16,   -7,    0],
       [  18,   10,    0],
       [  18,   18,    0],
       [  20,   16,    0],
       [  22,    1,    0],
       [  30,   12,    0],
       [  32,    5,    0],
       [  34,   18,    0],
       [  36,   10,    0],
       [  36,   23,    0],
       [  38,   -7,    0],
       [  38,    2,    0],
       [  38,    4,    0],
       [  38,   18,    0],
       [  50,   14,    0],
       [  52,   -2,    0],
       [  56,    5,    0],
       [  58,   28,    0],
       [  62,   11,    0],
       [  64,    6,    0],
       [  68,    6,    0],
       [  68,   16,    0],
       [  72,   -6,    0],
       [  78,   -5,    0],
       [  78,   19,    0],
       [  80,    2,    0],
       [  82,  -21,    0],
       [  82,   13,    0],
       [  82,   30,    0],
       [  84,    7,    0],
       [  86,   -6,    0],
       [  86,    5,    0],
       [  86,   23,    0],
       [  88,   27,    0],
       [  94,    7,    0],
       [  96,    8,    0],
       [  98,   13,    0],
       [ 102,   -3,    0],
       [ 106,   23,    0],
       [ 108,   18,    0],
       [ 112,   30,    0],
       [ 116,    6,    0],])

box = freud.box.Box(Lx=505, Ly=65, is2D=True)

voro = freud.locality.Voronoi()
voro.compute((box, centered_points));
nlist = voro.nlist
 
    
unwrapped_bonds = centered_points[nlist.query_point_indices] \
                    + box.wrap(centered_points[nlist.point_indices] \
                    - centered_points[nlist.query_point_indices])
real_connections = check_inside_box(unwrapped_bonds, box)

# just a check
slow_real_connections = box.contains(unwrapped_bonds)
assert(np.all(slow_real_connections==real_connections))

line_data = np.asarray(
    [[centered_points[i], centered_points[i] + box.wrap(centered_points[j] - centered_points[i])] for i, j in nlist[~real_connections]]
    )[:, :, :2]
line_collection = matplotlib.collections.LineCollection(line_data, alpha=0.5, colors ="red")

line_data_2 = np.asarray(
[[centered_points[i], centered_points[i] + box.wrap(centered_points[j] - centered_points[i])] for i, j in nlist[real_connections]]
)[:, :, :2]
line_collection_2 = matplotlib.collections.LineCollection(line_data_2, alpha=0.2, colors ="blue")

fig, ax=plt.subplots()
set_size(20,20)
ax = plt.gca()
voro.plot(ax=ax, cmap="RdBu")
ax.add_collection(line_collection)
ax.add_collection(line_collection_2)
plt.show()

[tommy-waltmann]

Upon closer inspection, it appears that the bonds that were thought to be missing are actually there, they just are in blue instead of red. This is a result of having a system which is short in the y-direction.

In the script posted above, all neighbor points which are separated by less than 1/2 box length in both the x and y directions have their bonds appear in blue, not red. Sometimes those points are voronoi neighbors of each other because one of their polytopes extends across the box. We would instinctive think those bonds should appear in red in the above plot, but they don't because of how short your system is in the y-direction.

Hope this clears everything up, let me know if there is any more confusion.

[Optimox]

@tommy-waltmann thanks for looking at it.

So is there a better way of detecting "artificial" bonds (or red bonds)?
If I understand correctly, if the two points are too close to each other in the y axis box.wrap(centered_points[j] - centered_points[i]) = centered_points[j] - centered_points[i] so they are not detected as artificial bonds ?

I understand that it is a difficult problem, but I can't lower the box size since my lower points are not on the same y level. That's why I started to ask if it was possible to use a convex hull instead of a box.

Do you see an effective way to discriminate between real and artificial bonds ?

[tommy-waltmann]

I think you have a few options going forward:

1. For your real system in 3D, it may be okay to use the method you have been using, depending on system parameters like box size, number of points, etc. The problem you are having is because of one short box dimension, if your 3D system doesn't have that, then you should be fine to proceed without changing anything.

2. I don't know what kind of system you are studying, but is it possible to just use a larger system at constant desnity of points? That would eliminate the problems you are having.

3. If you don't want to do either of those, then we have to do this the hard way. You will have to look at which voronoi polytopes cross the box, find the neighbors along the edges of the portion of the polytopes that cross the box, and make sure the bonds between those neighbors show up in red. The polytopes property of the Voronoi object will be helpful here. If the polytope of particle i crosses the box, then voro.polytopes[i] will have a vertex which is out of range of the normal box.

I would seriously recommend trying options 1 and 2 before going to option 3, but option 3 may still be doable if you play around with it long enough. If you tell me some more about what you are studying, I can be more helpful as far as recommending an option.

[Optimox]

@tommy-waltmann thanks for your help.

Unfortunately 1 and 2 are not possible for me.

However I ended up adding some fake points at top and bottom so that I just need to ignore connections from those fake points and end up with real connections only. Seems to be working ok. Thanks very much I'm closing this.