How to unfragment a "molecule"

Viewed 44

I have a "molecule" (which is really a graph consisting of 7 nodes) that is fragmented by periodic boundary conditions:

enter image description here

Essentially, the cell at the center is used to tile the entire space, and even though nodes 1 and 3 share an edge, they appear on the opposite sides of the cell at the center. I would like to unfragment this graph, to get:

enter image description here

Given

bonds =  [(0,1),(1,2),(1,3),(3,4),(3,5),(4,5),(5,6)]
coords = np.array([[4,5],[6,5.5],[6,4],[6.5,1],[1,2],[2,1],[3,1]])

and cell dimensions

cell = [7,6]

I try the following:

import networkx as nx
import numpy as np
import matplotlib.pyplot as plt

G = nx.Graph()
G.add_edges_from(bonds)

node = 2

processed = []
moved = []
Q = [node]
new_xyz = []
new_xyz.append(coords[node,:])
indices = [node]
while len(Q)>0:
    curr = Q[0]
    ref_pos = coords[curr,:]
    Q = Q[1:]
    processed.append(curr)
    for i in G.adj[curr]:
        if i not in processed:
            Q.append(i)
                if i not in moved:
                    print("moving ",i, "(ref: ",curr,")")
                    pos = coords[i,:]
                    dist = pos - ref_pos
                    dist = nearest_image(dist,cell)
                    pos2 = dist + ref_pos
                    new_xyz.append(pos2)
                    indices.append(i)
                    assert(curr in moved or curr==node)
                    moved.append(i)

new_xyz = np.array(new_xyz)

where

def nearest_image(d,cell):
    for i in range(2):
        hbox = cell[i] / 2

        while d[i] > hbox:
            d[i] -= cell[i]
        while d[i] < -hbox:
            d[i] += cell[i]
        assert(d[i] < hbox and d[i] > -hbox)

    return d

When I plot the results using:

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(5,3))
axes[0].scatter(coords[:,0],coords[:,1])
for i in range(7):
    axes[0].annotate(str(i), (coords[i,0],coords[i,1]))
axes[1].scatter(new_xyz[:,0],new_xyz[:,1])
for i in indices:
    axes[1].annotate(str(i), (new_xyz[i,0],new_xyz[i,1]))
plt.show()

I get:

enter image description here

What am I missing? I don't quite understand how 0 and 2 swap places either. The full script is accessible here.

1 Answers

If you detect that a vertex u is too far from a vertex v on the horizontal axis, then slide the whole figure by cell[0] / 2 horizontally. These two vertices are "too far" if they share an edge but are more than half the cell size apart.

Likewise for the vertical axis.

import networkx as nx
import matplotlib.pyplot as plt

bonds =  [(0,1),(1,2),(1,3),(3,4),(3,5),(4,5),(5,6)]
coords = [[4,5],[6,5.5],[6,4],[6.5,1],[1,2],[2,1],[3,1]]
cell = (7, 6)

# DRAW GRAPH WITH INITIAL POSITIONS
pos = dict(enumerate(coords))
G = nx.Graph()
G.add_edges_from(bonds)
nx.draw(G, pos=pos, with_labels=True)
plt.show()

# REPOSITION
dx = (cell[0] / 2) * any( 2 * abs(pos[u][0] - pos[v][0]) > cell[0] for u,v in bonds )
dy = (cell[1] / 2) * any( 2 * abs(pos[u][1] - pos[v][1]) > cell[1] for u,v in bonds )

new_pos = { v: ((x + dx) % cell[0], (y + dy) % cell[1]) for v,(x,y) in pos.items() }

# DRAW GRAPH WITH NEW POSITIONS
nx.draw(G, pos=new_pos, with_labels=True)
plt.show()

Initial positions: graph with initial position

After repositioning: graph after repositioning

Related