Home > Software design >  How to unfragment a "molecule"
How to unfragment a "molecule"

Time:07-28

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. graph with initial position

After repositioning: graph after repositioning

  • Related