Nearest Neighbors

One of the basic computations required for higher-level computations (such as the hexatic order parameter) is finding the nearest neighbors of a particle. This tutorial will show you how to compute the nearest neighbors and visualize that data.

The algorithm is straightforward:

for each particle i:
    for each particle j in neighbor_cells(i):
        r_ij = position[j] - position[i]
        r = sqrt(dot(r_ij, r_ij))
        l_r_array.append(r)
        l_n_array.append(j)
        # sort by distance
        sort(n_array, r_array)
        neighbor_array[i] = n_array[:k]

The data sets used in this example are a system of hard hexagons, simulated in the NVT thermodynamic ensemble in HOOMD-blue, for a dense fluid of hexagons at packing fraction \(\phi = 0.65\) and a solid at packing fractions \(\phi = 0.75\).

from bokeh.io import output_notebook
output_notebook()
from bokeh.models import Legend
from bokeh.plotting import figure, output_file, show
from bokeh.layouts import gridplot
import numpy as np
import time
from freud import parallel, box, locality
parallel.setNumThreads(4)
import util

# Create hexagon vertices
verts = util.make_polygon(sides=6, radius=0.6204)

# Define colors for our system
c_list = ["#30A2DA", "#FC4F30", "#E5AE38", "#6D904F", "#9757DB",
          "#188487", "#FF7F00", "#9A2C66", "#626DDA", "#8B8B8B"]
c_dict = dict()
c_dict[6] = c_list[0]
c_dict[5] = c_list[1]
c_dict[4] = c_list[2]
c_dict[3] = c_list[7]
c_dict[7] = c_list[4]

def render_plot(p, fbox):
    # Display box
    corners = util.box_2d_to_points(fbox)
    p.patches(xs=[corners[:-1, 0]], ys=[corners[:-1, 1]],
              fill_color=(0, 0, 0, 0), line_color="black", line_width=2)
    p.legend.location = 'bottom_center'
    p.legend.orientation = 'horizontal'
    util.default_bokeh(p)
    show(p)

Loading BokehJS ...

# Load the data
data_path = "data/phi065"
box_data = np.load("{}/box_data.npy".format(data_path))
pos_data = np.load("{}/pos_data.npy".format(data_path))
quat_data = np.load("{}/quat_data.npy".format(data_path))
n_frames = pos_data.shape[0]

Viewing our system

Before proceeding, we should probably view our system first. freud does not make any assumptions about your data and is not specifically designed for any one visualization package. Here we use bokeh to render our system. Bokeh is not appropriate for real-time interaction with your simulation data, nor is it appropriate for 3D data, but is perfectly fine for rendering individual simulation frames, so we will use it here

# Grab data from last frame
l_box = box_data[-1].tolist()
l_pos = pos_data[-1]
l_quat = quat_data[-1]
l_ang = 2*np.arctan2(np.copy(l_quat[:, 3]), np.copy(l_quat[:, 0]))

# Create box
fbox = box.Box.from_box(l_box)
side_length = max(fbox.Lx, fbox.Ly)
l_max = side_length / 2.0
l_max *= 1.1

# Take local vertices and rotate, translate into system coordinates
patches = util.local_to_global(verts, l_pos[:, :2], l_ang)

# Plot
p = figure(title="System Visualization",
           x_range=(-l_max, l_max),
           y_range=(-l_max, l_max))
p.patches(xs=patches[:, :, 0].tolist(), ys=patches[:, :, 1].tolist(),
          fill_color=(42, 126, 187), line_color="black", line_width=1.5, legend="Hexagons")
render_plot(p, fbox)

By eye, we can see regions where the hexagons appear close-packed, as well as regions where there are vacancies. We will be using the Nearest Neighbor object to investigate this in our system.

The Nearest Neighbor object

This module will give the indices of the \(k\) particles which are nearest to another particle. freud provides two different modes by which to compute the nearest neighbors, selected by the strict_cut variable:

• strict_cut=False (default): The value for rmax is expanded until every particle has \(k\) nearest neighbors

• strict_cut=True: the rmax value is not expanded, so that any “vacancies” in the number of neighbors found are filled with UINTMAX

strict_cut=False

First we show how to use the strict_cut=False mode to find the neighbors of a specific particle

# Create freud nearest neighbor object
n_neigh = 6
nn = locality.NearestNeighbors(rmax=1.5, n_neigh=n_neigh, strict_cut=False)

# Compute nearest neighbors
nn.compute(fbox, l_pos, l_pos)
# Get the NeighborList
n_list = nn.nlist
# Get the neighbors for particle 0
pidx = 0
n_idxs = n_list.index_j[np.where(n_list.index_i == pidx)[0]]

# Get position, orientation for the central particle
center_pos = l_pos[np.newaxis, pidx]
center_ang = l_ang[np.newaxis, pidx]

# Get the positions, orientations for the neighbor particles
neigh_pos = np.zeros(shape=(n_neigh, 3), dtype=np.float32)
neigh_ang = np.zeros(shape=(n_neigh), dtype=np.float32)
neigh_pos[:] = l_pos[n_idxs]
neigh_ang[:] = l_ang[n_idxs]

# Create array of transformed positions
c_patches = util.local_to_global(verts, center_pos[:, 0:2], center_ang)
n_patches = util.local_to_global(verts, neigh_pos[:, 0:2], neigh_ang)

# Create array of colors
center_color = np.array([c_list[0] for _ in range(center_pos.shape[0])])
neigh_color = np.array([c_list[-1] for _ in range(neigh_pos.shape[0])])

# Plot
p = figure(title="Nearest Neighbors Visualization",
           x_range=(-l_max, l_max), y_range=(-l_max, l_max))
p.patches(xs=n_patches[:, :, 0].tolist(), ys=n_patches[:, :, 1].tolist(),
          fill_color=neigh_color.tolist(), line_color="black", legend="Neighbors")
p.patches(xs=c_patches[:, :, 0].tolist(), ys=c_patches[:, :, 1].tolist(),
          fill_color=center_color.tolist(), line_color="black", legend="Center")
render_plot(p, fbox)

Notice that nearest neighbors properly handles periodic boundary conditions.

We do the same thing below, but for a particle/neighbors not spanning the box.

# Get the neighbors for particle 1000
pidx = 1000
segment_start = n_list.segments[pidx]
segment_end = segment_start + n_list.neighbor_counts[pidx]
n_idxs = n_list.index_j[segment_start:segment_end]

# Get position, orientation for the central particle
center_pos = l_pos[np.newaxis, pidx]
center_ang = l_ang[np.newaxis, pidx]

# Get positions, orientations for the neighbors and one non-neighbor
neigh_pos = l_pos[n_idxs]
neigh_ang = l_ang[n_idxs]
non_neigh_pos = neigh_pos[np.newaxis, -1]
non_neigh_ang = neigh_ang[np.newaxis, -1]

# Create array of transformed positions
c_patches = util.local_to_global(verts, center_pos[:, :2], center_ang)
n_patches = util.local_to_global(verts, neigh_pos[:, :2], neigh_ang)
non_n_patches = util.local_to_global(verts, non_neigh_pos[:, :2], non_neigh_ang)

# Create array of colors
center_color = np.array([c_list[0] for _ in range(center_pos.shape[0])])
neigh_color = np.array([c_list[-1] for _ in range(neigh_pos.shape[0])])

# Color the last particle differently
non_neigh_color = np.array([c_list[-2] for _ in range(non_neigh_pos.shape[0])])

# Plot
p = figure(title="Nearest Neighbors Visualization",
           x_range=(-l_max, l_max), y_range=(-l_max, l_max))
p.patches(xs=n_patches[:, :, 0].tolist(), ys=n_patches[:, :, 1].tolist(),
          fill_color=neigh_color.tolist(), line_color="black", legend="Neighbors")
p.patches(xs=non_n_patches[:, :, 0].tolist(), ys=non_n_patches[:, :, 1].tolist(),
          fill_color=non_neigh_color.tolist(), line_color="black", legend="Non-neighbor")
p.patches(xs=c_patches[:, :, 0].tolist(), ys=c_patches[:, :, 1].tolist(),
          fill_color=center_color.tolist(), line_color="black", legend="Center")
render_plot(p, fbox)

Notice that while freud found the 6 nearest neighbors, one of the particles isn’t really in a neighbor position (which we have colored purple). How do we go about finding particles with a deficit or surplus of neighbors?

strict_cut=True

Now for strict_cut=True. This mode allow you to find particles which have fewer than the specified number of particles. For this system, we’ll search for 8 neighbors, so that we can display particles with both a deficit and a surplus of neighbors.

# Create freud nearest neighbors object
n_neigh = 8
rmax = 1.7
nn = locality.NearestNeighbors(rmax=rmax, n_neigh=n_neigh, strict_cut=True)

# Compute nearest neighbors
nn.compute(fbox, l_pos, l_pos)
# Get the neighborlist
n_list = nn.nlist
# Get the number of particles
num_particles = nn.n_ref

p = figure(title="Nearest Neighbors visualization",
           x_range=(-l_max, l_max), y_range=(-l_max, l_max))
for k in np.unique(n_list.neighbor_counts):
    # Find particles with k neighbors
    c_idxs = np.where(n_list.neighbor_counts == k)[0]
    center_pos = l_pos[c_idxs]
    center_ang = l_ang[c_idxs]
    c_patches = util.local_to_global(verts, center_pos[:, 0:2], center_ang)
    center_color = np.array([c_dict[k] for _ in range(center_pos.shape[0])])
    p.patches(xs=c_patches[:, :, 0].tolist(), ys=c_patches[:, :, 1].tolist(),
              fill_color=center_color.tolist(), line_color="black", legend="k={}".format(k))
render_plot(p, fbox)

Visualize each set of k values independently

for k in np.unique(n_list.neighbor_counts):
    p = figure(title="Nearest Neighbors: k={}".format(k),
               x_range=(-l_max, l_max), y_range=(-l_max, l_max))
    # Find particles with k neighbors
    c_idxs = np.where(n_list.neighbor_counts == k)[0]
    center_pos = l_pos[c_idxs]
    center_ang = l_ang[c_idxs]
    c_patches = util.local_to_global(verts, center_pos[:, 0:2], center_ang)
    patches = util.local_to_global(verts, l_pos[:, 0:2], l_ang)
    center_color = np.array([c_dict[k] for _ in range(center_pos.shape[0])])
    p.patches(xs=patches[:, :, 0].tolist(), ys=patches[:, :, 1].tolist(),
              fill_color=(128, 128, 128, 0.1), line_color=(0, 0, 0, 0.1), legend="other")
    p.patches(xs=c_patches[:, :, 0].tolist(), ys=c_patches[:, :, 1].tolist(),
              fill_color=center_color.tolist(), line_color="black", legend="k={}".format(k))
    render_plot(p, fbox)

for k in np.unique(n_list.neighbor_counts):
    p = figure(title="Nearest Neighbors: k={}".format(k),
               x_range=(-l_max, l_max), y_range=(-l_max, l_max))
    # Find particles with k neighbors
    c_idxs = np.copy(np.where(n_list.neighbor_counts == k)[0])
    center_pos = l_pos[c_idxs]
    center_ang = l_ang[c_idxs]

    neigh_pos = np.zeros(shape=(k*len(c_idxs), 3), dtype=np.float32)
    neigh_ang = np.zeros(shape=(k*len(c_idxs)), dtype=np.float32)
    for i, pidx in enumerate(c_idxs):
        # Create a list of positions, angles to draw
        segment_start = n_list.segments[pidx]
        segment_end = segment_start + n_list.neighbor_counts[pidx]
        n_idxs = n_list.index_j[segment_start:segment_end]
        for j, nidx in enumerate(n_idxs):
            neigh_pos[k*i+j] = l_pos[nidx]
            neigh_ang[k*i+j] = l_ang[nidx]
    c_patches = util.local_to_global(verts, center_pos[:, 0:2], center_ang)
    n_patches = util.local_to_global(verts, neigh_pos[:, 0:2], neigh_ang)
    patches = util.local_to_global(verts, l_pos[:, 0:2], l_ang)
    center_color = np.array([c_dict[k] for _ in range(center_pos.shape[0])])
    neigh_color = np.array([c_list[-1] for _ in range(neigh_pos.shape[0])])
    p.patches(xs=patches[:, :, 0].tolist(), ys=patches[:, :, 1].tolist(),
              fill_color=(128, 128, 128, 0.1), line_color=(0, 0, 0, 0.1), legend="other")
    p.patches(xs=n_patches[:, :, 0].tolist(), ys=n_patches[:, :, 1].tolist(),
              fill_color=neigh_color.tolist(), line_color="black", legend="neighbors")
    p.patches(xs=c_patches[:, :, 0].tolist(), ys=c_patches[:, :, 1].tolist(),
              fill_color=center_color.tolist(), line_color="black", legend="k={}".format(k))
    render_plot(p, fbox)

Compare against a solid/crystalline system

Visualize in the same way, but with a solid system. Notice that almost all the hexagons have 6 nearest neighbors.

data_path = "data/phi075"
box_data = np.load("{}/box_data.npy".format(data_path))
pos_data = np.load("{}/pos_data.npy".format(data_path))
quat_data = np.load("{}/quat_data.npy".format(data_path))
n_frames = pos_data.shape[0]

# Grab data from last frame
l_box = box_data[-1].tolist()
l_pos = pos_data[-1]
l_quat = quat_data[-1]
l_ang = 2*np.arctan2(np.copy(l_quat[:, 3]), np.copy(l_quat[:, 0]))

# Create box
fbox = box.Box.from_box(l_box)
side_length = max(fbox.Lx, fbox.Ly)
l_max = side_length / 2.0
l_max *= 1.1

# Create freud nearest neighbor object
n_neigh = 8
rmax = 1.65
nn = locality.NearestNeighbors(rmax=rmax, n_neigh=n_neigh, strict_cut=True)

# Compute nearest neighbors
nn.compute(fbox, l_pos, l_pos)
# Get the neighborlist
n_list = nn.nlist

p = figure(title="Nearest Neighbors Visualization",
           x_range=(-l_max, l_max), y_range=(-l_max, l_max))
for k in np.unique(n_list.neighbor_counts):
    # find particles with k neighbors
    c_idxs = np.where(n_list.neighbor_counts == k)[0]
    center_pos = l_pos[c_idxs]
    center_ang = l_ang[c_idxs]
    c_patches = util.local_to_global(verts, center_pos[:, 0:2], center_ang)
    center_color = np.array([c_dict[k] for _ in range(center_pos.shape[0])])
    p.patches(xs=c_patches[:, :, 0].tolist(), ys=c_patches[:, :, 1].tolist(),
              fill_color=center_color.tolist(), line_color="black", legend="k={}".format(k))
render_plot(p, fbox)