Read standardized molecules as graphs

from pathlib import Path
import base64, gzip, hashlib, json, platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML
import networkx as nx
from IPython.display import display as _display
def display(value):
    if isinstance(value,pd.DataFrame):
        table=value.to_html(max_rows=12,escape=True)
        style="<style>.vibeit-readable-table{max-width:100%;overflow-x:auto}.vibeit-readable-table table{width:max-content;max-width:none;border-collapse:collapse}.vibeit-readable-table td,.vibeit-readable-table th{white-space:nowrap!important;word-break:normal!important;overflow-wrap:normal!important}</style>"
        return _display(HTML(style+'<div class="vibeit-readable-table">'+table+'</div>'))
    return _display(value)
RECIPE_ID='ch02-molecular-graphs'
OUTPUT_DIR=Path.cwd()/"results"/RECIPE_ID
plt.rcParams.update({"figure.dpi":120,"font.size":10,"axes.spines.top":False,"axes.spines.right":False,
                      "axes.prop_cycle":plt.cycler(color=["#18785f","#416aa6","#bc6541","#9974a5"])})
print("Python:",platform.python_version(),"; NumPy:",np.__version__,"; Pandas:",pd.__version__)
print("Working folder:",Path.cwd())
print("Core data mode: embedded snapshot; no network requests")

Python: 3.13.13 ; NumPy: 2.5.3 ; Pandas: 3.0.6
Working folder: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads
Core data mode: embedded snapshot; no network requests

def load_snapshot(recipe_id, expected_sha256, snapshot_key):
    """Read embedded data from this notebook, including a renamed copy."""
    for path in sorted(Path.cwd().glob("*.ipynb")):
        try:
            notebook = json.loads(path.read_text(encoding="utf-8"))
            metadata = notebook.get("metadata", {}).get("vibeit_cookbook", {})
            blob = metadata.get("snapshots", {}).get(snapshot_key)
            if metadata.get("recipe_id") != recipe_id or not blob:
                continue
            if blob.get("sha256") != expected_sha256:
                continue
            raw = gzip.decompress(base64.b64decode(blob["gzip_base64"], validate=True))
            if hashlib.sha256(raw).hexdigest() != expected_sha256:
                raise ValueError("Embedded snapshot checksum mismatch: " + snapshot_key)
            return json.loads(raw)
        except (OSError, json.JSONDecodeError):
            continue
    raise FileNotFoundError(
        "Save/open the complete .ipynb in the active working folder. "
        "The embedded snapshot could not be located for " + recipe_id)

SNAPSHOT_HASHES={'molecules': '780c5f98bb3b7261465f8330d4b31589d4224cdf21f916b4a1788f11bcc00aac'}
def snapshot(key):
    return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))

Embedded snapshots: ['molecules']

def parse_molfile(text):
    """Read V2000 atoms/bonds and formal-charge annotations from one SDF record."""
    lines=text.splitlines()
    if len(lines)<4 or 'V2000' not in lines[3]:raise ValueError('Expected one V2000 molfile/SDF record')
    atom_count=int(lines[3][:3]);bond_count=int(lines[3][3:6])
    if atom_count<=0 or len(lines)<4+atom_count+bond_count:raise ValueError('Truncated molfile')
    atoms=[];charge_codes={0:0,1:3,2:2,3:1,5:-1,6:-2,7:-3}
    for j,line in enumerate(lines[4:4+atom_count]):
        element=line[31:34].strip();position=[float(line[:10]),float(line[10:20]),float(line[20:30])]
        if not element or not np.isfinite(position).all():raise ValueError('Invalid atom record')
        charge_code=int(line[36:39].strip() or 0)
        if charge_code not in charge_codes:raise ValueError('Unsupported atom charge/radical code')
        atoms.append(dict(index=j,element=element,position=position,formal_charge=charge_codes[charge_code]))
    bonds=[];seen=set()
    for line in lines[4+atom_count:4+atom_count+bond_count]:
        a=int(line[:3])-1;b=int(line[3:6])-1;order=int(line[6:9]);pair=tuple(sorted([a,b]))
        if not 0<=a<atom_count or not 0<=b<atom_count or a==b or pair in seen or order not in [1,2,3,4]:
            raise ValueError('Invalid or duplicate bond')
        seen.add(pair);bonds.append(dict(a=a,b=b,order=order))
    for line in lines[4+atom_count+bond_count:]:
        if line.startswith('M  CHG'):
            fields=line.split();count=int(fields[2])
            for j in range(count):
                index=int(fields[3+2*j])-1;charge=int(fields[4+2*j])
                if not 0<=index<atom_count:raise ValueError('Invalid formal-charge atom index')
                atoms[index]['formal_charge']=charge
    return dict(atoms=atoms,bonds=bonds,title=lines[0],atom_count=atom_count,bond_count=bond_count)

print("Method ready: parse_molfile")

Method ready: parse_molfile

def molecular_graph(parsed):
    graph=nx.Graph()
    for atom in parsed['atoms']:graph.add_node(atom['index'],element=atom['element'],formal_charge=atom['formal_charge'])
    for bond in parsed['bonds']:graph.add_edge(bond['a'],bond['b'],reported_bond_type=bond['order'])
    components=nx.number_connected_components(graph)
    return graph,dict(atoms=graph.number_of_nodes(),bonds=graph.number_of_edges(),connected_components=components,cycle_rank=graph.number_of_edges()-graph.number_of_nodes()+components)

print("Method ready: molecular_graph")

Method ready: molecular_graph

molecules=snapshot("molecules")
properties=pd.DataFrame(molecules["properties"]).set_index("CID")
assert len(properties)==8 and properties.index.is_unique
display(properties[["name","MolecularFormula","MolecularWeight","TPSA","Charge"]])
parsed={int(cid):parse_molfile(molecules["sdf_2d"][str(cid)]) for cid in properties.index}
graphs={};summaries=[]
for cid,record in parsed.items():
 graph,summary=molecular_graph(record);graphs[cid]=graph
 summaries.append({"CID":cid,"name":properties.loc[cid,"name"],**summary,"formal_charge":sum(a["formal_charge"] for a in record["atoms"])})
summary=pd.DataFrame(summaries).set_index("CID");display(summary)

cid=241;record=parsed[cid];graph=graphs[cid]
positions={a["index"]:a["position"][:2] for a in record["atoms"]}
labels={a["index"]:a["element"]+str(a["index"]) for a in record["atoms"]}
fig,ax=plt.subplots(figsize=(7,6));nx.draw_networkx(graph,pos=positions,labels=labels,node_size=600,node_color=["#b4d4ed" if graph.nodes[j]["element"]=="C" else "#eeeeee" for j in graph],ax=ax,font_size=9)
nx.draw_networkx_edge_labels(graph,positions,edge_labels=nx.get_edge_attributes(graph,"reported_bond_type"),ax=ax,font_size=8)
ax.set_title("Benzene: reported 2D molecular graph");ax.axis("off");plt.tight_layout();plt.show()

fig,ax=plt.subplots(figsize=(8,4));ax.barh(summary.name,summary.cycle_rank)
ax.set(xlabel="Independent cycle rank",title="Topology from reported atom/bond records");plt.tight_layout();plt.show()

degrees=np.concatenate([np.array([degree for node,degree in graph.degree()]) for graph in graphs.values()])
fig,ax=plt.subplots(figsize=(7,4));ax.hist(degrees,bins=np.arange(-.5,6.6,1),edgecolor="white")
ax.set(xlabel="Reported atom-neighbor degree",ylabel="Atoms",title="Degree counts across declared structures");plt.tight_layout();plt.show()

for cid,record in parsed.items():
 assert graphs[cid].number_of_nodes()==record["atom_count"] and graphs[cid].number_of_edges()==record["bond_count"]
 assert summary.loc[cid,"formal_charge"]==properties.loc[cid,"Charge"]
assert summary.loc[241,"cycle_rank"]==1
OUTPUT_DIR.mkdir(parents=True,exist_ok=True);summary.to_csv(OUTPUT_DIR/"molecular-graph-summary.csv")
nx.write_graphml(graphs[241],OUTPUT_DIR/"benzene.graphml")
print("Header, graph and charge checks passed; exports:",OUTPUT_DIR)

Header, graph and charge checks passed; exports: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads/results/ch02-molecular-graphs

RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
    import requests
    try:
        response=requests.get('https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/cid/2519/property/MolecularFormula/JSON',timeout=20)
        response.raise_for_status()
        print("Source:",response.url)
        print(response.text[:700])
    except requests.RequestException as error:
        print("Online extension unavailable:",type(error).__name__,str(error)[:160])
else:
    print("Online extension skipped; embedded core data unchanged")

Online extension skipped; embedded core data unchanged