Implementing Common Neighbor Analysis as a custom methodΒΆ

Researchers commonly wish to implement their own custom analysis methods for particle simulations. Here, we show an example of how to write Common Neighbor Analysis (Honeycutt and Andersen, J. Phys. Chem. 91, 4950) as a custom method using freud and the NetworkX package.

NetworkX can be installed with pip install networkx.

First, we generate random points and determine which points share neighbors.

[1]:
from collections import defaultdict

import freud
import numpy as np
[2]:
# Use a face-centered cubic (fcc) system
box, points = freud.data.UnitCell.fcc().generate_system(4)
aq = freud.AABBQuery(box, points)
nl = aq.query(points, {"num_neighbors": 12, "exclude_ii": True}).toNeighborList()
[3]:
# Get all sets of common neighbors.
common_neighbors = defaultdict(list)
for i, p in enumerate(points):
    for j in nl.point_indices[nl.query_point_indices == i]:
        for k in nl.point_indices[nl.query_point_indices == j]:
            if i != k:
                common_neighbors[(i, k)].append(j)

Next, we use NetworkX to build graphs of common neighbors and compute the Common Neighbor Analysis signatures.

[4]:
from collections import Counter

import networkx as nx

diagrams = defaultdict(list)
particle_counts = defaultdict(Counter)

for (a, b), neighbors in common_neighbors.items():
    # Build up the graph of connections between the
    # common neighbors of a and b.
    g = nx.Graph()
    for i in neighbors:
        for j in set(nl.point_indices[nl.query_point_indices == i]).intersection(
            neighbors
        ):
            g.add_edge(i, j)

    # Define the identifiers for a CNA diagram:
    # The first integer is 1 if the particles are bonded, otherwise 2
    # The second integer is the number of shared neighbors
    # The third integer is the number of bonds among shared neighbors
    # The fourth integer is an index, just to ensure uniqueness of diagrams
    diagram_type = 2 - int(b in nl.point_indices[nl.query_point_indices == a])
    key = (diagram_type, len(neighbors), g.number_of_edges())
    # If we've seen any neighborhood graphs with this signature,
    # we explicitly check if the two graphs are identical to
    # determine whether to save this one. Otherwise, we add
    # the new graph immediately.
    if key in diagrams:
        isomorphs = [nx.is_isomorphic(g, h) for h in diagrams[key]]
        if any(isomorphs):
            idx = isomorphs.index(True)
        else:
            diagrams[key].append(g)
            idx = diagrams[key].index(g)
    else:
        diagrams[key].append(g)
        idx = diagrams[key].index(g)
    cna_signature = key + (idx,)
    particle_counts[a].update([cna_signature])

Looking at the counts of common neighbor signatures, we see that the first particle of the fcc structure has 12 bonds with signature \((1, 4, 2, 0)\) as we expect.

[5]:
particle_counts[0]
[5]:
Counter({(1, 4, 2, 0): 12,
         (2, 4, 4, 0): 6,
         (2, 1, 0, 0): 12,
         (2, 2, 1, 0): 24})