Biochar Structure Visualization Demo
This notebook demonstrates generating and visualizing biochar molecular structures.
Requirements (in addition to ``biochar``):
conda install -c conda-forge py3dmol jupyter matplotlib pandas
[ ]:
import py3Dmol
import numpy as np
import matplotlib.pyplot as plt
from rdkit.Chem import Draw, MolToMolBlock
from rdkit.Chem.Draw import rdMolDraw2D
from IPython.display import display, SVG
from biochar import generate_biochar, generate_surface, generate_biochar_series
Helper: 3D viewer
A small utility that builds a py3Dmol view from an RDKit mol + coords array. Using the coords array (rather than the mol’s internal conformer) ensures the viewer always reflects the final post-processed positions.
[ ]:
def coords_to_xyz(mol, coords):
"""Build an XYZ-format string from an RDKit mol and a coords array (Å)."""
lines = [str(mol.GetNumAtoms()), ""]
for atom, (x, y, z) in zip(mol.GetAtoms(), coords):
lines.append(f"{atom.GetSymbol()} {x:.4f} {y:.4f} {z:.4f}")
return "\n".join(lines)
def view3d(mol, coords, width=700, height=450, style="stick"):
"""Display an interactive 3D structure."""
v = py3Dmol.view(width=width, height=height)
v.addModel(coords_to_xyz(mol, coords), "xyz")
if style == "stick":
v.setStyle({"stick": {"radius": 0.12},
"sphere": {"scale": 0.25}})
elif style == "sphere":
v.setStyle({"sphere": {"colorscheme": "Jmol", "scale": 0.4}})
# Colour oxygen red so functional groups stand out
v.setStyle({"elem": "O"}, {"sphere": {"color": "red", "scale": 0.45},
"stick": {"color": "red", "radius": 0.15}})
v.zoomTo()
return v
1 Single molecule — 2D structure
[ ]:
mol, coords, gro, top, itp = generate_biochar(
target_num_carbons=50,
H_C_ratio=0.5,
O_C_ratio=0.1,
functional_groups={"phenolic": 3, "ether": 1},
output_directory="output/demo",
basename="bc50",
seed=42,
)
print(f"{mol.GetNumAtoms()} atoms — {mol.GetNumHeavyAtoms()} heavy")
# 2D depiction (kekulized for clarity)
drawer = rdMolDraw2D.MolDraw2DSVG(600, 380)
drawer.drawOptions().addAtomIndices = False
drawer.DrawMolecule(mol)
drawer.FinishDrawing()
display(SVG(drawer.GetDrawingText()))
2 Single molecule — interactive 3D
[ ]:
view3d(mol, coords).show()
3 Functional group comparison
Side-by-side 2D depictions of the same carbon skeleton with different oxygen functional groups.
[ ]:
fg_configs = [
("Phenolic (–OH)", {"phenolic": 4}),
("Carboxyl (–COOH)", {"carboxyl": 2}),
("Ether (–O–)", {"ether": 3}),
("Mixed", {"phenolic": 2, "carboxyl": 1, "ether": 1}),
]
mols, labels = [], []
for label, fg in fg_configs:
m, *_ = generate_biochar(
target_num_carbons=38,
functional_groups=fg,
output_directory="output/demo/fg",
basename=label.split()[0].lower(),
seed=7,
)
mols.append(m)
labels.append(label)
img = Draw.MolsToGridImage(mols, molsPerRow=2, subImgSize=(500, 320),
legends=labels)
display(img)
4 Pentagon ring defects
Topological disorder increases with defect_fraction; the structures become non-planar at higher fractions.
[ ]:
defect_mols, defect_coords, defect_labels = [], [], []
for frac in [0.0, 0.10, 0.20]:
m, c, *_ = generate_biochar(
target_num_carbons=60,
defect_fraction=frac,
output_directory="output/demo/defects",
basename=f"defect_{int(frac*100):02d}pct",
seed=99,
)
defect_mols.append(m)
defect_coords.append(c)
defect_labels.append(f"defect_fraction={frac:.2f}")
display(Draw.MolsToGridImage(defect_mols, molsPerRow=3, subImgSize=(420, 280),
legends=defect_labels))
[ ]:
# 3D view of the 20% defect structure — notice the out-of-plane curvature
view3d(defect_mols[2], defect_coords[2]).show()
5 Slit-pore surface
Two parallel graphene-like sheets separated by a user-defined pore diameter.
[ ]:
sheets, gro_surf, top_surf, itp_surf = generate_surface(
target_num_carbons=40,
functional_groups={"phenolic": 2, "ether": 1},
pore_diameter=10.0,
output_directory="output/demo/surface",
basename="slit_10A",
seed=42,
)
print(f"Generated {len(sheets)} sheets")
for i, s in enumerate(sheets):
zmin, zmax = s.coords[:, 2].min(), s.coords[:, 2].max()
print(f" Sheet {i+1}: {s.mol.GetNumAtoms()} atoms, "
f"z ∈ [{zmin:.1f}, {zmax:.1f}] Å")
[ ]:
# Visualise both sheets in one viewer, coloured differently
SHEET_COLOURS = ["#4C8EDA", "#E07B39"] # blue / orange
v = py3Dmol.view(width=750, height=500)
for i, sheet in enumerate(sheets):
v.addModel(coords_to_xyz(sheet.mol, sheet.coords), "xyz")
v.setStyle(
{"model": i},
{"stick": {"color": SHEET_COLOURS[i], "radius": 0.12},
"sphere": {"color": SHEET_COLOURS[i], "scale": 0.25}},
)
# Highlight oxygens in red regardless of sheet colour
v.setStyle(
{"model": i, "elem": "O"},
{"sphere": {"color": "red", "scale": 0.45},
"stick": {"color": "red", "radius": 0.15}},
)
v.zoomTo()
v.show()
6 Temperature series — van Krevelen diagram
A van Krevelen plot (H/C vs O/C) is the standard way to characterise biochar across pyrolysis temperatures. Lower temperature → upper-right; higher temperature → lower-left.
[ ]:
from biochar import BiocharGenerator, GeneratorConfig
temp_configs = [
{"label": "300–400 °C", "target_num_carbons": 70, "H_C_ratio": 0.70, "O_C_ratio": 0.22, "seed": 1},
{"label": "400–500 °C", "target_num_carbons": 80, "H_C_ratio": 0.60, "O_C_ratio": 0.16, "seed": 2},
{"label": "500–600 °C", "target_num_carbons": 90, "H_C_ratio": 0.50, "O_C_ratio": 0.11, "seed": 3},
{"label": "600–700 °C", "target_num_carbons": 100, "H_C_ratio": 0.38, "O_C_ratio": 0.06, "seed": 4},
{"label": "700–800 °C", "target_num_carbons": 110, "H_C_ratio": 0.28, "O_C_ratio": 0.03, "seed": 5},
]
results = []
for cfg in temp_configs:
gen = BiocharGenerator(GeneratorConfig(
target_num_carbons=cfg["target_num_carbons"],
H_C_ratio=cfg["H_C_ratio"],
O_C_ratio=cfg["O_C_ratio"],
seed=cfg["seed"],
))
_, _, comp = gen.generate()
results.append({
"label": cfg["label"],
"H_C_target": cfg["H_C_ratio"],
"O_C_target": cfg["O_C_ratio"],
"H_C_actual": comp.H_C_ratio,
"O_C_actual": comp.O_C_ratio,
"n_C": comp.num_carbons,
})
[ ]:
fig, ax = plt.subplots(figsize=(7, 5))
cmap = plt.cm.plasma
colours = [cmap(i / (len(results) - 1)) for i in range(len(results))]
for r, c in zip(results, colours):
ax.scatter(r["O_C_actual"], r["H_C_actual"], s=120, color=c,
edgecolors="k", linewidths=0.6, zorder=3)
ax.annotate(r["label"], (r["O_C_actual"], r["H_C_actual"]),
textcoords="offset points", xytext=(8, 4), fontsize=9)
# Target values as ×
for r, c in zip(results, colours):
ax.scatter(r["O_C_target"], r["H_C_target"], marker="x", s=60,
color=c, linewidths=1.5, zorder=2)
ax.set_xlabel("O/C atomic ratio", fontsize=12)
ax.set_ylabel("H/C atomic ratio", fontsize=12)
ax.set_title("Van Krevelen diagram — simulated biochar series", fontsize=13)
ax.legend(
handles=[
plt.scatter([], [], s=80, c="gray", edgecolors="k", label="Actual"),
plt.scatter([], [], marker="x", s=60, c="gray", label="Target"),
],
fontsize=9,
)
ax.grid(True, alpha=0.3)
ax.set_xlim(-0.02, 0.30)
ax.set_ylim(0.0, 0.90)
plt.tight_layout()
plt.show()
7 Asymmetric slit-pore surface
Different chemistry on each wall — e.g. one oxidised wall and one clean wall.
[ ]:
asym_sheets, *_ = generate_surface(
pore_diameter=12.0,
sheet_overrides=[
{"target_num_carbons": 38, "functional_groups": {"phenolic": 4, "carboxyl": 1}, "seed": 10},
{"target_num_carbons": 38, "functional_groups": {"ether": 2}, "seed": 11},
],
output_directory="output/demo/asymmetric",
basename="asym_12A",
)
# 2D depictions of each wall
wall_mols = [s.mol for s in asym_sheets]
wall_labels = ["Wall 1: phenolic + carboxyl", "Wall 2: ether"]
display(Draw.MolsToGridImage(wall_mols, molsPerRow=2, subImgSize=(500, 320),
legends=wall_labels))
[ ]:
v2 = py3Dmol.view(width=750, height=500)
for i, sheet in enumerate(asym_sheets):
v2.addModel(coords_to_xyz(sheet.mol, sheet.coords), "xyz")
v2.setStyle(
{"model": i},
{"stick": {"color": SHEET_COLOURS[i], "radius": 0.12},
"sphere": {"color": SHEET_COLOURS[i], "scale": 0.25}},
)
v2.setStyle(
{"model": i, "elem": "O"},
{"sphere": {"color": "red", "scale": 0.45},
"stick": {"color": "red", "radius": 0.15}},
)
v2.zoomTo()
v2.show()