Inspect a calculated caffeine conformer

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 py3Dmol
from scipy.spatial.distance import cdist
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='ch03-caffeine-geometry'
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

molecules=snapshot("molecules");sdf=molecules["caffeine_3d_sdf"];record=parse_molfile(sdf)
positions=np.array([a["position"] for a in record["atoms"]]);elements=np.array([a["element"] for a in record["atoms"]])
heavy=elements!="H";print("Atoms:",len(positions),"; heavy atoms:",int(heavy.sum()),"; coordinate unit: Å; computed conformer")
viewer=py3Dmol.view(width="100%",height=420)
viewer.addModel(sdf,"sdf");viewer.setStyle({"stick":{"radius":.15},"sphere":{"scale":.25}});viewer.setBackgroundColor("white");viewer.zoomTo();viewer.zoom(.7)
display(HTML(viewer._make_html()))

Atoms:

 24 ; heavy atoms: 14 ; coordinate unit: Å; computed conformer

colors={"C":"gray","H":"lightgray","N":"royalblue","O":"red"}
fig=plt.figure(figsize=(7,6));ax=fig.add_subplot(projection="3d")
for element in sorted(set(elements)):
 mask=elements==element;ax.scatter(*positions[mask].T,color=colors.get(element,"green"),s=35 if element!="H" else 12,label=element)
for bond in record["bonds"]:
 pair=positions[[bond["a"],bond["b"]]];ax.plot(*pair.T,color="black",lw=1)
ax.set(xlabel="x (Å)",ylabel="y (Å)",zlabel="z (Å)",title="Calculated caffeine conformer; source bonds");ax.legend();plt.tight_layout();plt.show()

bond_lengths=np.array([np.linalg.norm(positions[b["a"]]-positions[b["b"]]) for b in record["bonds"]])
fig,ax=plt.subplots(figsize=(7,4));ax.hist(bond_lengths,bins=15,edgecolor="white")
ax.set(xlabel="Calculated source-bond length (Å)",ylabel="Bonds",title="Geometric lengths, not interaction energies");plt.tight_layout();plt.show()

heavy_positions=positions[heavy];distance=cdist(heavy_positions,heavy_positions)
fig,ax=plt.subplots(figsize=(6,5));im=ax.imshow(distance,cmap="viridis")
ax.set(xlabel="Heavy-atom order",ylabel="Heavy-atom order",title="Calculated heavy-atom distances");fig.colorbar(im,ax=ax,label="Distance (Å)");plt.tight_layout();plt.show()

assert np.allclose(distance,distance.T) and np.allclose(np.diag(distance),0)
assert np.isclose(distance[0,1],np.sqrt(np.sum((heavy_positions[0]-heavy_positions[1])**2)))
assert np.all(bond_lengths>0) and heavy.sum()==14
OUTPUT_DIR.mkdir(parents=True,exist_ok=True);(OUTPUT_DIR/"caffeine-computed-conformer.sdf").write_text(sdf)
pd.DataFrame(distance).to_csv(OUTPUT_DIR/"heavy-atom-distances-Angstrom.csv",index=False)
pd.DataFrame({"a":[b["a"] for b in record["bonds"]],"b":[b["b"] for b in record["bonds"]],"reported_type":[b["order"] for b in record["bonds"]],"length_Angstrom":bond_lengths}).to_csv(OUTPUT_DIR/"reported-bond-lengths.csv",index=False)
print("Geometry checks passed; exports:",OUTPUT_DIR)

Geometry checks passed; exports: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads/results/ch03-caffeine-geometry

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