"""Write report/network.html, a clickable map of every species and step the explorer found, from data/ only.
python reproduce/network_page.py # run from the Finding's folder
Needs numpy. The page is one file; its 3D viewer (3Dmol.js) loads from cdnjs.cloudflare.com.
"""
import csv
import json
from pathlib import Path
import numpy as np
HERE = Path(__file__).resolve().parent.parent
PARTIAL = (1.25, 1.6) # bonds being made or broken at a TS: 1.25-1.6x the sum of covalent radii
RCOV = {"H": 0.31, "C": 0.76, "O": 0.66} # covalent radii, Angstrom (Cordero et al. 2008)
def rows(name):
with open(HERE / "data" / f"{name}.csv", newline="", encoding="utf-8") as fh:
return list(csv.DictReader(fh))
def frames(fp):
lines, out, i = fp.read_text().split("\n"), {}, 0
while i < len(lines) and lines[i].strip():
n = int(lines[i])
out[lines[i + 1].strip()] = "\n".join(lines[i:i + 2 + n]) + "\n"
i += 2 + n
return out
def stretched_bonds(xyz_text):
"""Atom pairs of a TS whose bond is being made or broken: stretched, not sharing a bonded neighbour,
not H-H."""
body = xyz_text.split("\n")[2:]
sym = [ln.split()[0] for ln in body if ln.strip()]
x = np.array([[float(v) for v in ln.split()[1:4]] for ln in body if ln.strip()])
rc = np.array([RCOV[e] for e in sym])
r = np.linalg.norm(x[:, None] - x[None], axis=-1) / (rc[:, None] + rc[None])
n = len(sym)
tight = r < PARTIAL[0]
np.fill_diagonal(tight, False)
return [[j, i] for i in range(n) for j in range(i) if PARTIAL[0] <= r[i, j] < PARTIAL[1]
and not (tight[i] & tight[j]).any() and (sym[i] != "H" or sym[j] != "H")]
def main():
species, steps = rows("species"), rows("steps")
xyz = frames(HERE / "data" / "network.xyz")
num = lambda v: None if v == "" else float(v)
g = {0: 0.0}
changed = True # free energies of the species the kinetics reach, from the open steps' barriers
while changed:
changed = False
for st in steps:
if st["in_reach"] != "1":
continue
a, b = int(st["from_species"]), int(st["to_species"])
ts = float(st["ts_dG_kcal"])
for k, bar in ((a, st["dG_barrier_fwd_kcal"]), (b, st["dG_barrier_rev_kcal"])):
if k not in g:
g[k] = ts - float(bar)
changed = True
sp_out = []
for r in species:
k = int(r["species"])
sp_out.append({"id": k, "smiles": r["smiles"], "stereo": r["stereo"], "e": float(r["rel_energy_kcal"]),
"g": round(g[k], 3) if k in g else None, "frac": float(r["final_fraction"]),
"xyz": xyz[f"species {k}"]})
st_out = []
for r in steps:
i = int(r["step"])
st_out.append({"id": i, "a": int(r["from_species"]), "b": int(r["to_species"]), "kind": r["kind"],
"ets": float(r["ts_rel_energy_kcal"]), "bf": float(r["barrier_fwd_kcal"]),
"br": float(r["barrier_rev_kcal"]), "gts": num(r["ts_dG_kcal"]),
"gbf": num(r["dG_barrier_fwd_kcal"]), "gbr": num(r["dG_barrier_rev_kcal"]),
"open": r["in_reach"] == "1", "xyz": xyz[f"step {i} transition state"],
"made": stretched_bonds(xyz[f"step {i} transition state"])})
data = json.dumps({"species": sp_out, "steps": st_out}, separators=(",", ":"))
page = TEMPLATE.replace("/*DATA*/", data.replace("
Cyclopentadiene + Methyl Acrylate Network
endoexoin reachout of reach
Click a species (bar) or a transition state (‡) for its structure and energies.