get the interface along x-axis for 2D active matter system. Avoid internal bubbles. Very cool.

interface.py

import re
import freud
import numpy as np
import networkx as nx
from io import StringIO
from scipy.spatial.distance import squareform, pdist


def should_join(p1, p2):
    return 0 in pdist(np.concatenate((p1, p2))[:, None])


def join_pairs(pairs, copy=True):
    """
    Args:
        pairs (list): a list of tuples
        copy (bool): very difficult
    Example:
        >>> pairs = [(2, 3), (3, 5), (2, 6), (8, 9), (9, 10)]
        >>> join_pairs(pairs)
        [(2, 3, 5, 6), (8, 9, 10)]
    """
    pairs_joined = pairs.copy() if copy else pairs
    for p1 in pairs_joined:
        for p2 in pairs_joined:
            if (p1 is not p2) and should_join(p1, p2):
                pairs_joined.append(tuple(set(p1 + p2)))
                pairs_joined.remove(p1)
                pairs_joined.remove(p2)
                pairs_joined = join_pairs(pairs_joined, copy=False)
                return pairs_joined  # join once per recursion
    return pairs_joined

def load_str_from_file(f, rows):
    """
    load data from a file handler to numpy array
    Args:
        f (TextIOWrapper): a file handler obtained from `open`
        rows (int): the number of rows to load into the numpy array
    """
    data_str = ""
    for _ in range(rows):
        data_str += f.readline()
    io_stream = StringIO(data_str)
    return np.loadtxt(io_stream)


def parse_dump(filename):
    """
    Args:
        filename (str): path to the .dump file
    Return:
        dict: the configuration and metadata in different frames
    """
    f = open(filename, 'r')
    frames = []
    for line in f:
        if re.match(r'ITEM: TIMESTEP\n', line):
            time_step = int(re.search('\d+', f.readline()).group(0))
            f.readline()  # jump through one line
            particle_num = int(re.search('\d+', f.readline()).group(0))
            f.readline()  # jump through one line
            box = load_str_from_file(f, 3)
            col_names = re.split('\s+', f.readline())[1:]
            configuration = load_str_from_file(f, particle_num)
            frames.append({
                "conf": configuration,
                "col_names": col_names,
                "box": box,
                "particle_num": particle_num,
                "time_step": time_step,
            })
    f.close()
    return frames


def pdist_pbc(positions, box):
    """
    Get the pair-wise distances of particles in a priodic boundary box
    Args:
        positiosn (:obj:`numpy.ndarray`): coordinate of particles, shape (N, dim)
        box (:obj:`numpy.ndarray`): the length of the priodic bondary, shape (dim,)
    Return:
        :obj:`numpy.ndarray`: the pairwise distance, shape ( (N * N - N) / 2, ),
            use :obj:`scipy.spatial.distance.squareform` to recover the matrix form.
    """
    n, dim = positions.shape
    result = np.zeros(int((n * n - n) / 2), dtype=np.float64)
    for d in range(dim):
        dist_1d = pdist(positions[:, d][:, np.newaxis])
        dist_1d[dist_1d > box[d] / 2] -= box[d]
        result += np.power(dist_1d, 2)
    return np.sqrt(result)


def get_interface_from_voro(points, box, threshold):
    """
    Use the voronoi volume to get the interface particles
    Args:
        points (np.ndarray): shape (n, 3), the position of z axis
            should all be zero for freud
        box (np.ndarray): the box size of x and y axis, shape (2,)
        threshold (double): if the volume is above threshold value
        the particle is the interface
    Return:
        np.ndarray: the points corresponding to the interface, may contain bubbles
    """
    voro = freud.locality.Voronoi()
    box_obj = freud.box.Box.from_box(box)
    box_obj.periodic_x = True
    box_obj.periodic_y = True
    voro.compute((box_obj, points))
    vv = voro.volumes
    interface = points[vv > threshold, :2]
    return interface


def get_percolating_cycles(
    points, box, neighbour_distance, percolate_cutoff=0.9, connection_cutoff=4
):
    """
    Generate a graph based on distance threshold and find cycles in the graph
        that is percolating in the X-axis

    Args:
        points (np.ndarray): shape (n, 2), the position of the interface
            that may contain bubboles
        box (np.ndarray): the pbc box size of x and y axis, shape (2,)
        neighbour_distance (float): an edge is if two particles have distance
            smaller than this value
        percolate_cutoff (float): a percolate circle shoud span this_value * box_length
        connection_cutoff (double): a percolate circle should be continues, the minimum
            distance between two particles should be smaller than this_value * nn_distance
    """
    distances = squareform(pdist_pbc(points, box))  # adjacancy matrix
    adjacency = distances < neighbour_distance
    graph = nx.from_numpy_matrix(adjacency)

    cycles = nx.cycle_basis(graph)
    cycles_percolate = []

    num_cutoff = (box[0] * percolate_cutoff) / neighbour_distance

    for c in cycles:
        if len(c) > num_cutoff:
            indices = np.array(c)
            x_vals = points[indices, 0]
            x_sorted = np.sort(x_vals)
            is_percolate = (x_vals.max() - x_vals.min()) > box[0] * percolate_cutoff
            no_gap = (x_sorted[1:] - x_sorted[:-1]).max() < neighbour_distance * connection_cutoff
            if is_percolate and no_gap:
                cycles_percolate.append(indices)
    for i, cp in enumerate(cycles_percolate):
        neighbours = np.concatenate([
            np.fromiter(graph.neighbors(node), dtype=int) for node in cp
        ])
        cycles_percolate[i] = np.concatenate((cp, neighbours), axis=0)
    return cycles_percolate


def merge_overlap_interface_indices(interface_indices):
    """
    Merge the interface indices if they have a common element/index
    """
    n = len(interface_indices)
    if n <= 1:
        return interface_indices
    else:
        should_merge = []
        for i in range(n):
            for j in range(i + 1, n):
                if len(np.intersect1d(
                    interface_indices[i], interface_indices[j]
                )) > 0:
                    should_merge.append((i, j))
        if len(should_merge) == 0:
            return interface_indices
        else:
            to_merge = join_pairs(should_merge)
            result = [interface_indices[x] for x in range(n) if x not in np.concatenate(to_merge)]
            for indices in to_merge:
                merged = np.unique(np.concatenate([interface_indices[x] for x in indices]))
                result.append(merged)
            return result


def get_interface_refined(
    points, box, neighbour_distance, percolate_cutoff=0.9, connection_cutoff=4
):
    """
    Remove deceptive interface points that were actually internal bubbles
        by finding circles in a Graph representation of the data
    Args:
        points (np.ndarray): shape (n, 2), the position of the interface
            that may contain bubboles
        box (np.ndarray): the pbc box size of x and y axis, shape (2,)
    Return:
        np.ndarray: the points corresponding to the true interfaces
    """
    interface_indices = get_percolating_cycles(
        points, box, neighbour_distance, percolate_cutoff, connection_cutoff
    )
    merged_indices = merge_overlap_interface_indices(interface_indices)
    return [points[mi] for mi in merged_indices]

def get_dnn_mean(points, box):
    """
    Get the average nearest neighbour distances in a PBC box
    Args:
        points (np.ndarray): shape (n, 2), the position of the interface
            that may contain bubboles
        box (np.ndarray): the pbc box size of x and y axis, shape (2,)
    Return:
        np.ndarray: the points corresponding to the true interfaces
    """
    distances = squareform(pdist_pbc(points, box))
    np.fill_diagonal(distances, np.inf)
    return np.min(distances, axis=0).mean()


if __name__ == "__main__":
    import matplotlib.pyplot as plt

    dump_file = '2000.lammpstrj'

    frames = parse_dump(dump_file)
    for fn, f in enumerate(frames[1:]):
        xyz = f['conf'][:, 2:5]
        plt.scatter(*xyz.T[:2], alpha=0.1, color='k')
        box = f['box'][:2, 1]
        psi_6 = f['conf'][:, 15]
        interface = get_interface_from_voro(xyz, box, 1.0)  # 2D system, Z=0
        plt.scatter(*interface.T, color='k', alpha=0.5, facecolor='none', label='Voro Interface')
        dnn = get_dnn_mean(interface, box)
        interfaces = get_interface_refined(interface, box, dnn * 2.0)
        for interface in interfaces:
            plt.scatter(*interface.T, marker='+', label='Refined Interface')
        plt.xlabel("X", fontsize=16)
        plt.ylabel("Y", fontsize=16)
        plt.legend(fontsize=16)
        plt.xlim(0, box[0])
        plt.ylim(0, box[1])
        plt.gcf().set_size_inches(box[0]/10, box[1]/10)
        plt.savefig(f'frame-{fn:04}.png')
        plt.close()

This is what I got with the above code.

This is what I got from the 2nd Version. The different interfaces were correctly identified if the interfaces were not too bad.