Problem making graph structures of COFs

[aahil1412]

Hi, I am trying to get a graph representation of covalent organic frameworks (COF). However, I faced a problem while creating the plot. So, COFs have periodic structures like crystals, and I want to incorporate them into the graph structure.

I am attaching the code I used below:

from pymatgen.core import Structure
from pymatgen.io.cif import CifParser
import networkx as nx
import plotly.graph_objects as go
# Step 1: Read the CIF file using the updated parse_structures method
cif_file = 'COF-5.cif'  # Replace with your actual CIF file path
parser = CifParser(cif_file)
structure = parser.parse_structures(primitive=False)[0]  # Get the first structure from the CIF file
# Step 2: Create a graph from the structure, with adjusted cutoffs and periodicity handling
def structure_to_graph(structure, cutoff=None, use_periodicity=False):
    """
    Converts a pymatgen Structure object to a NetworkX graph.    
    Parameters:
    - structure: pymatgen Structure object
    - cutoff: dict, cutoff distances for different atom pairs (in angstroms)
    - use_periodicity: bool, whether to consider periodic boundary conditions   
    Returns:
    - G: NetworkX Graph object
    """
    G = nx.Graph()   
    # Add nodes with features
    for i, site in enumerate(structure.sites):
        atom_type = site.specie.symbol
        atomic_number = site.specie.number
        coords = site.coords  # Cartesian coordinates
        G.add_node(i, 
                   atom=atom_type,
                   atomic_number=atomic_number,
                   coords=coords)
    # Add edges based on bonding criteria with specific distance cutoffs
    for i, site in enumerate(structure.sites):
        neighbors = structure.get_neighbors(site, r=2.0, include_image=use_periodicity)  # Adjust r and periodicity
        for neighbor in neighbors:
            j = neighbor.index
            if i < j:  # Avoid duplicate edges in an undirected graph             
                # Define cutoff distances for specific atom pairs (C-H, N-H, etc.)
                atom_pair = tuple(sorted([site.specie.symbol, neighbor.specie.symbol]))
                bond_cutoff = cutoff.get(atom_pair, 1.8)  # Default cutoff if not specified
                # Only add edges if the bond length is within the atom-specific cutoff
                if neighbor.nn_distance < bond_cutoff:
                    G.add_edge(i, j, distance=neighbor.nn_distance)
    
    return G
# Define cutoff distances for specific bonds (in angstroms)
cutoff_distances = {
    ('C', 'H'): 1.1,
    ('C', 'C'): 1.5,
    ('C', 'N'): 1.4,
    ('N', 'H'): 1.1,
    ('O', 'H'): 1.0,
}
# Convert the structure to a graph with bond-specific cutoffs and without periodicity
graph = structure_to_graph(structure, cutoff=cutoff_distances, use_periodicity=False)
# Step 3: Assign colors to atoms based on type
atom_color_map = {
    'C': 'blue',
    'O': 'red',
    'H': 'green',
    'N': 'purple',
    'Ni': 'orange',
    'S': 'yellow'
    # Add other atom types if needed
}
# Extract node positions directly from real coordinates
node_x = []
node_y = []
node_colors = []  # Store the colors for each node
for node in graph.nodes():
    coords = graph.nodes[node]['coords']  # Real 3D coordinates
    node_x.append(coords[0])  # Use X and Y for the 2D plot
    node_y.append(coords[1])
    # Get the atom type and assign the corresponding color
    atom_type = graph.nodes[node]['atom']
    node_colors.append(atom_color_map.get(atom_type, 'gray'))  # Default to gray if atom type is not in the map
# Create edge data based on real atomic positions
edge_x = []
edge_y = []
edge_text = []  # To store interatomic distances
for edge in graph.edges(data=True):
    x0, y0 = graph.nodes[edge[0]]['coords'][0], graph.nodes[edge[0]]['coords'][1]
    x1, y1 = graph.nodes[edge[1]]['coords'][0], graph.nodes[edge[1]]['coords'][1]
    edge_x.append(x0)
    edge_x.append(x1)
    edge_x.append(None)  # Separating each segment
    edge_y.append(y0)
    edge_y.append(y1)
    edge_y.append(None)   
    # Extract and format the bond distance
    distance = edge[2]['distance']
    edge_text.append(f"Distance: {distance:.2f} Å")
# Create edge trace with hover text showing distances
edge_trace = go.Scatter(
    x=edge_x, y=edge_y,
    line=dict(width=1, color='#888'),
    hoverinfo='text',
    mode='lines',
    text=edge_text  # Add hover information for distances
)
# Create node trace with different colors for atom types
node_trace = go.Scatter(
    x=node_x, y=node_y,
    mode='markers',
    hoverinfo='text',
    marker=dict(
        color=node_colors,  # Use the colors based on atom type
        size=10,
        colorscale='YlGnBu',  # Optional, can remove if you want only atom-based colors
        showscale=False  # Turn off color bar since we're using discrete colors
    ))
# Add node labels and connections as hover text
node_text = []
for node, adjacencies in enumerate(graph.adjacency()):
    atom_type = graph.nodes[node]['atom']
    node_text.append(f'Node {node}, Atom: {atom_type}, Connections: {len(adjacencies[1])}')
node_trace.text = node_text
# Customize layout for better spacing
fig = go.Figure(data=[edge_trace, node_trace],
                layout=go.Layout(
                    title='Network Graph with Bond-Specific Cutoffs and No Periodicity',
                    titlefont_size=16,
                    showlegend=False,
                    hovermode='closest',
                    margin=dict(b=0, l=0, r=0, t=40),
                    xaxis=dict(showgrid=False, zeroline=False),
                    yaxis=dict(showgrid=False, zeroline=False))
                )
# Show the interactive plot
fig.show()

This is the output I got:

However, the actual structure should look like this:

I tried using this code for a Metal-organic framework (which has a more defined crystal system than COFs), and it worked perfectly, but in this case, the periodic conditions mismatched.

I also tried doing it for another COf structure, but that led to the same result, too.

I am not sure what I need to change in the code. I am also attaching the cif file below.

Thank you
COF-5.zip

[DanielYang59]

Hi @aahil1412 I believe you are just one step from there!

You would need to check the distance between two sites (i and j) instead of the nn_distance (to nearest neighbor).

# Only add edges if the bond length is within the atom-specific cutoff
- if neighbor.nn_distance < bond_cutoff:
+ if structure.get_distance(i, j, jimage=0) < bond_cutoff:  # use jimage=0 to ignore PBC
    G.add_edge(i, j, distance=neighbor.nn_distance)

Then the plot becomes: