Back to blog

Drawing Reactions: Highlighting and Mapping Reactions with RDKit

rdkitcheminformaticsvisualizationpublication
JR
Joseph Rheinhardt, PhD
28 min read
Drawing Reactions: Highlighting and Mapping Reactions with RDKit

Highlighting and Mapping Reactions

This post extends the molecule-highlighting ideas from Part 3a: Highlighting and Annotating Molecules to reaction schemes. We will draw a DCC-mediated amide coupling, apply abbreviations to keep each panel readable, and then highlight the reaction center atoms in each component.

The example reaction

Our example is the DCC-mediated coupling of Boc-Trp-OH with H-Tyr-OMe to give Boc-Trp-Tyr-OMe; this is the same dipeptide used as the example molecule in Part 3a.

The reaction was sourced from De Novo Chem's Saguaro Chem reaction search engine. Searching for an acid-amine coupling of Boc-protected tryptophan with a tyrosine methyl ester returned this DCC/HOBt/NMM procedure in DCM with a reported 77% yield.

Saguaro Chem Reaction Details panel for the Boc-Trp-OH + H-Tyr-OMe coupling

Saguaro Chem makes it easy to copy the entire reaction SMILES or to copy the reactant and product SMILES from their molecule detail windows.

CC(C)(C)OC(=O)N[C@@H](Cc1c[nH]c2ccccc12)C(=O)O   (Boc-Trp-OH, R1)
+ COC(=O)[C@@H](N)Cc1ccc(O)cc1                     (H-Tyr-OMe, R2)
→ COC(=O)[C@H](Cc1ccc(O)cc1)NC(=O)[C@@H](Cc1c[nH]c2ccccc12)NC(=O)OC(C)(C)C

Reagents: DCC, HOBt, N-methylmorpholine; solvent: DCM; yield: 77%.

Setup

The helper functions used below are included in the Appendix at the end of this post. Copy them into a single script before running the examples.

from rdkit import Chem
from rdkit.Chem import AllChem

RXN_SMILES = (
    "CC(C)(C)OC(=O)N[C@@H](Cc1c[nH]c2ccccc12)C(=O)O"
    ".COC(=O)[C@@H](N)Cc1ccc(O)cc1"
    ">>COC(=O)[C@H](Cc1ccc(O)cc1)NC(=O)[C@@H](Cc1c[nH]c2ccccc12)NC(=O)OC(C)(C)C"
)

Baseline reaction scheme

A single draw_reaction call composes reactant and product panels around a labeled arrow:

draw_reaction(
    RXN_SMILES,
    "reaction_baseline.svg",
    agent_text="DCC, HOBt, N-methylmorpholine",
    condition_text="DCM",
)

Baseline reaction scheme

Abbreviations on reaction panels

Both reactants and the product contain fragments that can be condensed. Reactant 1 carries the NHBoc carbamate; Reactant 2 carries the 4-hydroxyphenyl (tyrosine) ring.

There is a labeling convention to observe: on the reactant side the phenol is a free hydroxyl group acting as a nucleophile, so the label reads PhOH (phenyl first, hydroxyl at the end). In the product the group is a substituent on the peptide backbone, so the convention from Part 3a applies: HOPh (hydroxyl closest to the bond). This means we need two separate abbreviation sets and must pass them independently to each side of the reaction:

reactant_abbrevs = parse_abbreviations("BocHN *NC(=O)OC(C)(C)C\nPhOH *c1ccc(O)cc1")
product_abbrevs  = parse_abbreviations("BocHN *NC(=O)OC(C)(C)C\nHOPh *c1ccc(O)cc1")

draw_reaction(
    RXN_SMILES,
    "reaction_abbreviated.svg",
    agent_text="DCC, HOBt, N-methylmorpholine",
    condition_text="DCM",
    reactant_abbreviations=reactant_abbrevs,
    product_abbreviations=product_abbrevs,
)

Abbreviated reaction scheme

Three implementation details

Getting abbreviations to work correctly on reaction panels required three implementation choices.

1. Coverage limit. RDKit's CondenseMolAbbreviations silently skips any abbreviation that would cover more than maxCoverage of the molecule's heavy atoms (default 40%). The tyrosine reactant R2 has only 14 heavy atoms; the PhOH/HOPh ring covers 8 of them (57%), which exceeds the default. The reaction helper introduces a dedicated function that always uses max_coverage=1.0, which is appropriate for the small reactant fragments that reactions typically contain:

def _apply_panel_abbreviations(mol: Chem.Mol, abbreviations) -> Chem.Mol:
    """Apply abbreviations to a reaction panel molecule.

    Uses ``max_coverage=1.0`` so that fragments such as PhOH/HOPh are not
    silently skipped on small reactant molecules where the ring would otherwise
    exceed the default 40% coverage limit.
    """
    return apply_abbreviations(mol, abbreviations, max_coverage=1.0)

2. SMILES round-trip. When a mol object has been abbreviated by RDKit, calling Chem.MolToSmiles() on it emits * wildcard atoms for each abbreviation attachment point. If that SMILES string is then passed into draw_highlighted_molecule, RDKit parses it fresh and renders the wildcards visibly in the image. The fix is to pass the already-abbreviated mol object directly to the drawing backend rather than converting it back to SMILES:

def _panel_svg(mol, highlight_config=None, abbreviations=None):
    if abbreviations is not None:
        mol = _apply_panel_abbreviations(mol, abbreviations)
    if highlight_config:
        return _draw_mol_highlighted_bytes(mol, highlight_config, fmt="svg").decode("utf-8")
    return draw_mol_bytes(mol, size=PANEL_SIZE, fmt="svg").decode("utf-8")

_draw_mol_highlighted_bytes is a new internal helper that accepts a mol object directly, builds the RDKit color maps, and calls PrepareAndDrawMolecule without any SMILES conversion.

3. Per-side abbreviations. A single abbreviations argument applied globally cannot produce different labels on the two sides of the arrow. draw_reaction was extended with reactant_abbreviations and product_abbreviations parameters that override the global abbreviations for their respective sides:

r_abbrevs = reactant_abbreviations if reactant_abbreviations is not None else abbreviations
p_abbrevs = product_abbreviations  if product_abbreviations  is not None else abbreviations

The molecule-level helpers are unchanged.

Highlighting reaction centers

To draw attention to the bond-forming event we highlight the reacting atoms in each component. The carboxylic acid of Reactant 1 (atoms 12, 13, 14 after abbreviation), the amine nitrogen of Reactant 2 (atom 5), and the new amide nitrogen and carbonyl in the product (atoms 7, 8, 9) are colored independently:

r1_atom_colors = {12: "#e74c3c", 13: "#e74c3c", 14: "#e74c3c"}
r2_atom_colors = {5: "#3498db"}
p_atom_colors  = {7: "#3498db", 8: "#e74c3c", 9: "#e74c3c"}

draw_reaction(
    RXN_SMILES,
    "reaction_highlighted_atoms.svg",
    agent_text="DCC, HOBt, N-methylmorpholine",
    condition_text="DCM",
    reactant_abbreviations=reactant_abbrevs,
    product_abbreviations=product_abbrevs,
    reactant_highlight_configs=[
        {"atom_colors": r1_atom_colors},
        {"atom_colors": r2_atom_colors},
    ],
    product_highlight_configs=[
        {"atom_colors": p_atom_colors},
    ],
)

Reaction with atom highlights on the reacting centers

A complementary view highlights only the newly formed C-N bond in the product, leaving the reactants unannotated:

p_bond_colors = {(7, 8): "#2ecc71"}

draw_reaction(
    RXN_SMILES,
    "reaction_highlighted_bond.svg",
    agent_text="DCC, HOBt, N-methylmorpholine",
    condition_text="DCM",
    reactant_abbreviations=reactant_abbrevs,
    product_abbreviations=product_abbrevs,
    product_highlight_configs=[
        {"bond_colors": p_bond_colors},
    ],
)

Reaction with the newly formed amide bond highlighted

Note that the atom indices used here are those of the abbreviated molecules, not the original SMILES. Running each component through draw_molecule with style={"addAtomIndices": True} before choosing indices is a reliable way to confirm them.

Atom-map numbers

A lightweight way to annotate which atoms participate in bond-forming and bond-breaking events is to attach small numeric labels — map numbers — directly to the atoms of interest. RDKit exposes this through the atomNote property, which draw_reaction accepts per panel via the atom_maps key in reactant_highlight_configs and product_highlight_configs (a dict mapping atom index to map number). Here we use map 1 for the carbonyl carbon, map 2 for the carbonyl oxygen, and map 3 for the amine nitrogen:

r1_atom_maps = {12: 1, 13: 2}   # carboxylic acid C=O
r2_atom_maps = {5: 3}            # free amine N
p_atom_maps  = {8: 1, 9: 2, 7: 3}  # new amide C=O and N

draw_reaction(
    RXN_SMILES,
    "reaction_atom_mapped.svg",
    agent_text=AGENTS,
    condition_text=CONDITIONS,
    reactant_abbreviations=reactant_abbrevs,
    product_abbreviations=product_abbrevs,
    reactant_highlight_configs=[
        {"atom_maps": r1_atom_maps},
        {"atom_maps": r2_atom_maps},
    ],
    product_highlight_configs=[
        {"atom_maps": p_atom_maps},
    ],
    style={"annotationFontScale": 0.7},
)

Full reaction with atom map numbers on the reaction centers

It is worth being clear about what these labels are and are not. They are manual annotations chosen to visually connect corresponding atoms across the reaction. They are not a true atom map in the cheminformatics sense. A genuine atom map requires solving the atom-mapping problem: determining, for every heavy atom in the product, exactly which atom in the reactants it came from, accounting for all possible bond-breaking and bond-forming pathways. That is a non-trivial combinatorial problem. De Novo Chem has built Agave Chem, a dedicated atom-mapping tool, to solve it properly. If you need rigorous atom maps for reaction analysis, synthesis planning, or machine learning, check out our open source page.

Appendix: complete helper code

Copy the following code into a single script before running the examples above.

# ---------- abbreviations helpers ----------
"""Chemical abbreviation helpers for RDKit drawings.

RDKit provides a built-in abbreviation engine in ``rdkit.Chem.rdAbbreviations``.
This module wraps it with a small set of common abbreviations and helpers so the
rest of the drawing toolkit can condense fragments before rendering without
knowing the details of the RDKit API.
"""
from rdkit import Chem
from rdkit.Chem import rdAbbreviations

# Common abbreviations used in medicinal-chemistry figures.  Each line is:
#   label SMARTS
# The ``*`` atom is the attachment point to the rest of the molecule.
DEFAULT_ABBREVIATIONS = """
Ac *C(=O)C
Allyl *CC=C
Bn *Cc1ccccc1
Boc *C(=O)OC(C)(C)C
Bz *C(=O)c1ccccc1
Cbz *C(=O)OCc1ccccc1
CF3 *C(F)(F)F
CHO *C=O
CO2Et *C(=O)OCC
COOH *C(=O)O
Cp *c1cccc1
Et *CC
Fmoc *C(=O)OCC1c2ccccc2-c3ccccc31
iPr *C(C)C
Me *C
Ms *S(=O)(=O)C
NMe *NC
NO2 *N(=O)=O
OMe *OC
OEt *OCC
Ph *c1ccccc1
tBu *C(C)(C)C
Tf *S(=O)(=O)C(F)(F)F
Ts *S(=O)(=O)c1ccc(C)cc1
"""


def parse_abbreviations(text: str) -> list:
    """Parse a text block of ``label SMARTS`` lines into RDKit abbreviations.

    Blank lines are ignored.  The returned object is a list of
    ``AbbreviationDefinition`` instances that can be passed to
    ``apply_abbreviations`` or ``CondenseMolAbbreviations``.
    """
    return rdAbbreviations.ParseAbbreviations(text)


def apply_abbreviations(
    mol: Chem.Mol,
    abbrev_defs,
    max_coverage: float = 0.4,
) -> Chem.Mol:
    """Return a copy of ``mol`` with matched fragments condensed to labels.

    Parameters
    ----------
    mol
        RDKit molecule to condense.
    abbrev_defs
        Abbreviation definitions (from ``parse_abbreviations`` or
        ``rdAbbreviations.GetDefaultAbbreviations``).
    max_coverage
        Maximum fraction of the molecule that a single abbreviation may cover
        before RDKit skips it.  Increase to 1.0 for very small molecules where
        the abbreviation is most of the structure.
    """
    mol_copy = Chem.Mol(mol)
    return rdAbbreviations.CondenseMolAbbreviations(
        mol_copy, abbrev_defs, maxCoverage=max_coverage
    )

# ---------- moldraw helpers ----------
"""Shared molecule drawing helpers, carried over from Post 2."""
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem.Draw import rdMolDraw2D

DEFAULT_STYLE = dict(
    fixedBondLength=30,
    bondLineWidth=2,
    minFontSize=14,
    maxFontSize=14,
    padding=0.05,
)


def _make_drawer(size, fmt):
    return rdMolDraw2D.MolDraw2DSVG(*size) if fmt == "svg" else rdMolDraw2D.MolDraw2DCairo(*size)


def draw_mol_bytes(
    mol,
    size=(300, 300),
    fmt: str = "png",
    style: dict = None,
    abbreviations=None,
) -> bytes:
    """Render an existing RDKit molecule object and return raw image bytes.

    This does **not** recompute 2D coordinates, so any alignment applied to
    ``mol`` before calling this function is preserved in the final image.

    If ``abbreviations`` is supplied, the molecule is copied and condensed
    before rendering so matched fragments are drawn as compact labels.
    """
    if abbreviations is not None:

        mol = apply_abbreviations(mol, abbreviations)

    drawer = _make_drawer(size, fmt)
    opts = drawer.drawOptions()
    for key, value in {**DEFAULT_STYLE, **(style or {})}.items():
        setattr(opts, key, value)
    opts.useBWAtomPalette()

    rdMolDraw2D.PrepareAndDrawMolecule(drawer, mol)
    drawer.FinishDrawing()
    data = drawer.GetDrawingText()
    return data.encode("utf-8") if fmt == "svg" else data


def draw_molecule_bytes(
    smiles: str,
    size=(300, 300),
    fmt: str = "png",
    style: dict = None,
    abbreviations=None,
) -> bytes:
    """Render a single molecule from SMILES and return the raw image bytes (PNG or SVG)."""
    mol = Chem.MolFromSmiles(smiles)
    AllChem.Compute2DCoords(mol)
    return draw_mol_bytes(mol, size=size, fmt=fmt, style=style, abbreviations=abbreviations)


def draw_molecule(
    smiles: str,
    path: str,
    size=(300, 300),
    fmt: str = None,
    style: dict = None,
    abbreviations=None,
) -> None:
    """Render a single molecule directly to a file.

    If ``fmt`` is not given, it is inferred from the file extension:
    ``.svg`` produces SVG, anything else produces PNG.

    If ``abbreviations`` is supplied, matched fragments are drawn as compact
    labels instead of full structures.
    """
    fmt = fmt or ("svg" if path.lower().endswith(".svg") else "png")
    data = draw_molecule_bytes(
        smiles, size=size, fmt=fmt, style=style, abbreviations=abbreviations
    )
    mode = "w" if fmt == "svg" else "wb"
    with open(path, mode) as f:
        f.write(data if fmt != "svg" else data.decode("utf-8"))

# ---------- highlight helpers ----------
"""Add highlighting, atom-mapping numbers, legends, and annotations to RDKit drawings.

This module extends the drawing helpers from ``the molecule renderer`` with visual emphasis:
atom-level highlighting, bond-level highlighting, clean reaction atom-map numbers,
and publication-style legends and annotations.  Both PNG and SVG output are
supported.
"""
import html
import io
import re
from typing import Iterable

from PIL import Image, ImageDraw, ImageFont
from rdkit import Chem, Geometry
from rdkit.Chem import AllChem
from rdkit.Chem.Draw import rdMolDraw2D, PrepareMolForDrawing


LEGEND_LINE_HEIGHT = 22
LEGEND_SWATCH_SIZE = 14
LEGEND_MARGIN = 10
FONT_SIZE = 14
ANNOTATION_HEIGHT = 22


def _infer_fmt(path: str, fmt: str | None = None) -> str:
    """Return 'svg' for ``.svg`` paths, otherwise 'png'."""
    if fmt:
        return fmt.lower()
    return "svg" if path.lower().endswith(".svg") else "png"


def _load_font(size: int):
    try:
        return ImageFont.truetype("DejaVuSans.ttf", size)
    except OSError:
        return ImageFont.load_default()


def _hex_to_rgb(color: str) -> tuple[float, float, float]:
    """Convert a hex color string to an RDKit (r, g, b) tuple in [0, 1]."""
    color = str(color).lstrip("#")
    if len(color) == 3:
        color = "".join(c * 2 for c in color)
    if len(color) != 6:
        raise ValueError(f"Invalid hex color: {color!r}")
    r, g, b = int(color[0:2], 16), int(color[2:4], 16), int(color[4:6], 16)
    return (r / 255.0, g / 255.0, b / 255.0)


def _to_rdkit_color_map(color_map: dict, mol: Chem.Mol) -> dict:
    """Convert user color keys (atom idx or atom pair) into an RDKit color map.

    Atom pairs are mapped to the bond index that connects them.  If no bond
    exists between the atoms, the entry is skipped.
    """
    rdkit_map: dict[int, tuple[float, float, float]] = {}
    for key, color in color_map.items():
        rdkit_color = _hex_to_rgb(color)
        if isinstance(key, int):
            rdkit_map[key] = rdkit_color
        elif isinstance(key, (tuple, list)) and len(key) == 2:
            bond = mol.GetBondBetweenAtoms(int(key[0]), int(key[1]))
            if bond is not None:
                rdkit_map[bond.GetIdx()] = rdkit_color
    return rdkit_map


def _bond_indices_from_atom_pairs(mol: Chem.Mol, pairs: Iterable[tuple[int, int]]) -> list[int]:
    """Return the RDKit bond indices for the given (atom1, atom2) pairs."""
    indices = []
    for a1, a2 in pairs:
        bond = mol.GetBondBetweenAtoms(int(a1), int(a2))
        if bond is not None:
            indices.append(bond.GetIdx())
    return indices


def _set_atom_map_notes(mol: Chem.Mol, atom_maps: dict[int, int]) -> None:
    """Attach RDKit ``atomNote`` properties so mapped atoms get numbered labels."""
    for atom_idx, map_num in atom_maps.items():
        atom = mol.GetAtomWithIdx(atom_idx)
        atom.SetProp("atomNote", str(map_num))


def _draw_shaded_rings(
    drawer,
    mol: Chem.Mol,
    ring_shades: list[tuple[list[int], str]],
    highlight_atoms: list[int] | None,
    bond_indices: list[int],
    highlight_atom_colors: dict[int, tuple[float, float, float]],
    highlight_bond_colors: dict[int, tuple[float, float, float]],
) -> None:
    """Draw ring interiors as shaded polygons underneath the molecule.

    This is a two-pass trick: RDKit must be drawn once to finalize atom
    coordinates and canvas scale, then the canvas is cleared, the polygons
    are filled, and the molecule is redrawn with ``clearBackground=False``
    so the shading shows through.
    """
    mol = PrepareMolForDrawing(mol)

    # First pass: set coordinates and scale.
    drawer.DrawMolecule(
        mol,
        highlightAtoms=highlight_atoms or [],
        highlightBonds=bond_indices,
        highlightAtomColors=highlight_atom_colors,
        highlightBondColors=highlight_bond_colors,
    )

    # Fill the interior of each requested ring.
    conf = mol.GetConformer()
    for ring_atoms, color in ring_shades:
        rdkit_color = _hex_to_rgb(color)
        pts = [
            Geometry.Point2D(conf.GetAtomPosition(aidx).x, conf.GetAtomPosition(aidx).y)
            for aidx in ring_atoms
        ]
        drawer.SetFillPolys(True)
        drawer.SetColour(rdkit_color)
        drawer.DrawPolygon(pts)

    # Second pass: redraw molecule on top of the shaded polygons.
    opts = drawer.drawOptions()
    opts.clearBackground = False
    drawer.DrawMolecule(
        mol,
        highlightAtoms=highlight_atoms or [],
        highlightBonds=bond_indices,
        highlightAtomColors=highlight_atom_colors,
        highlightBondColors=highlight_bond_colors,
    )


def _parse_rdkit_svg(svg_text: str) -> tuple[str, str]:
    """Extract the viewBox (or width/height) and inner content from an RDKit molecule SVG."""
    viewbox_match = re.search(r"viewBox=(?:'|\")([^'\"]+)(?:'|\")", svg_text)
    if viewbox_match:
        viewbox = viewbox_match.group(1)
    else:
        width_match = re.search(r"width=(?:'|\")(\d+)(?:'|\")", svg_text)
        height_match = re.search(r"height=(?:'|\")(\d+)(?:'|\")", svg_text)
        w = int(width_match.group(1)) if width_match else 100
        h = int(height_match.group(1)) if height_match else 100
        viewbox = f"0 0 {w} {h}"

    tag_start = svg_text.find("<svg")
    start = svg_text.find(">", tag_start)
    end = svg_text.rfind("</svg>")
    content = svg_text[start + 1 : end].strip() if start != -1 and end != -1 else ""
    return viewbox, content


def _legend_height(legend_items: list[tuple[str, str]]) -> int:
    """Return the vertical space required for the legend, or 0 if empty."""
    if not legend_items:
        return 0
    return LEGEND_MARGIN * 2 + len(legend_items) * LEGEND_LINE_HEIGHT


def _annotation_height(annotation: str | None) -> int:
    """Return the vertical space required for an annotation, or 0 if absent."""
    return ANNOTATION_HEIGHT + 8 if annotation else 0


def _draw_png_legend(
    draw: ImageDraw.ImageDraw,
    legend_items: list[tuple[str, str]],
    width: int,
    y_start: int,
) -> None:
    """Draw color swatches and labels on a PIL canvas."""
    font = _load_font(FONT_SIZE)
    x = LEGEND_MARGIN
    y = y_start + LEGEND_MARGIN
    for color, label in legend_items:
        draw.rectangle(
            [(x, y), (x + LEGEND_SWATCH_SIZE, y + LEGEND_SWATCH_SIZE)],
            fill=color,
            outline="black",
        )
        draw.text((x + LEGEND_SWATCH_SIZE + 6, y), label, font=font, fill="black")
        y += LEGEND_LINE_HEIGHT


def _draw_png_annotation(
    draw: ImageDraw.ImageDraw,
    annotation: str,
    width: int,
    y: int,
) -> None:
    """Draw a centered annotation line on a PIL canvas."""
    font = _load_font(FONT_SIZE)
    tw = draw.textlength(annotation, font=font)
    draw.text(((width - tw) / 2, y), annotation, font=font, fill="black")


def _svg_legend_items(
    legend_items: list[tuple[str, str]],
    x_offset: int,
    y_start: int,
) -> str:
    """Return SVG ``<rect>`` and ``<text>`` elements for the legend."""
    parts = []
    y = y_start + LEGEND_MARGIN
    for color, label in legend_items:
        parts.append(
            f'<rect x="{x_offset}" y="{y}" width="{LEGEND_SWATCH_SIZE}" '
            f'height="{LEGEND_SWATCH_SIZE}" fill="{color}" stroke="black" stroke-width="1"/>'
        )
        parts.append(
            f'<text x="{x_offset + LEGEND_SWATCH_SIZE + 6}" y="{y + LEGEND_SWATCH_SIZE - 1}" '
            f'font-size="{FONT_SIZE}" font-family="sans-serif">{html.escape(label)}</text>'
        )
        y += LEGEND_LINE_HEIGHT
    return "\n".join(parts)


def _svg_annotation(annotation: str, width: int, y: int) -> str:
    """Return an SVG ``<text>`` element centered below the panel."""
    return (
        f'<text x="{width / 2}" y="{y + FONT_SIZE}" text-anchor="middle" '
        f'font-size="{FONT_SIZE}" font-family="sans-serif">{html.escape(annotation)}</text>'
    )


def draw_highlighted_molecule(
    smiles: str,
    path: str,
    size: tuple[int, int] = (400, 400),
    highlight_atoms: list[int] | None = None,
    highlight_bonds: list[tuple[int, int]] | None = None,
    atom_colors: dict[int, str] | None = None,
    bond_colors: dict[tuple[int, int], str] | None = None,
    ring_shades: list[tuple[list[int], str]] | None = None,
    atom_maps: dict[int, int] | None = None,
    legend: dict[str, str] | None = None,
    annotation: str | None = None,
    fmt: str | None = None,
    abbreviations=None,
    max_coverage: float = 0.4,
    style: dict | None = None,
) -> None:
    """Draw a single molecule with optional highlights, atom maps, legend, and annotation.

    Parameters
    ----------
    smiles
        SMILES string for the molecule.
    path
        Output file path.  Format is inferred from the extension unless ``fmt``
        is supplied.
    size
        Width and height of the molecule panel (legend/annotation add extra
        canvas height below the panel).
    highlight_atoms
        Atom indices to highlight.
    highlight_bonds
        (atom1, atom2) pairs defining bonds to highlight.
    atom_colors
        Mapping from atom index to hex color for per-atom highlight colors.
    bond_colors
        Mapping from (atom1, atom2) pair to hex color for per-bond highlight colors.
    ring_shades
        Optional list of ``(atom_indices, color)`` tuples that shade the interior
        of each specified ring. This uses a two-pass drawing trick because
        RDKit does not natively support filled ring shading.
    atom_maps
        Mapping from atom index to map number, drawn as a small atom note.
    legend
        Mapping from hex color to label, shown as swatches below the panel.
    annotation
        Optional caption text centered below the panel.
    fmt
        Output format override (``"png"`` or ``"svg"``).
    abbreviations
        Optional abbreviation definitions.  If supplied, matched fragments are
        condensed into labels before rendering.
    max_coverage
        Maximum fraction of the molecule that a single abbreviation may cover
        before RDKit skips it.  Default is 0.4; raise to 1.0 for small molecules.
    style
        Optional drawing options that override ``DEFAULT_STYLE`` (e.g.
        ``{"highlightRadius": 0.3}``).
    """
    fmt = _infer_fmt(path, fmt)
    mol = Chem.MolFromSmiles(smiles)
    AllChem.Compute2DCoords(mol)

    highlight_atom_colors = _to_rdkit_color_map(atom_colors or {}, mol)
    highlight_bond_colors = _to_rdkit_color_map(bond_colors or {}, mol)
    bond_indices = _bond_indices_from_atom_pairs(mol, highlight_bonds or [])

    # If the caller supplied color maps but no explicit highlight lists,
    # assume every colored atom/bond should be highlighted.
    if highlight_atoms is None and atom_colors:
        highlight_atoms = list(atom_colors.keys())
    if highlight_bonds is None and bond_colors:
        bond_indices = _bond_indices_from_atom_pairs(mol, list(bond_colors.keys()))

    panel_w, panel_h = size
    legend_items = list(legend.items()) if legend else []
    legend_h = _legend_height(legend_items)
    annot_h = _annotation_height(annotation)
    total_w = panel_w
    total_h = panel_h + legend_h + annot_h

    if abbreviations is not None:

        mol = apply_abbreviations(mol, abbreviations, max_coverage=max_coverage)

    if atom_maps:
        _set_atom_map_notes(mol, atom_maps)

    drawer = _make_drawer((panel_w, panel_h), fmt)
    opts = drawer.drawOptions()
    for key, value in {**DEFAULT_STYLE, **(style or {})}.items():
        setattr(opts, key, value)
    opts.useBWAtomPalette()

    if ring_shades:
        _draw_shaded_rings(
            drawer,
            mol,
            ring_shades,
            highlight_atoms,
            bond_indices,
            highlight_atom_colors,
            highlight_bond_colors,
        )
    else:
        rdMolDraw2D.PrepareAndDrawMolecule(
            drawer,
            mol,
            highlightAtoms=highlight_atoms or [],
            highlightBonds=bond_indices,
            highlightAtomColors=highlight_atom_colors,
            highlightBondColors=highlight_bond_colors,
        )
    drawer.FinishDrawing()

    if fmt == "svg":
        viewbox, content = _parse_rdkit_svg(drawer.GetDrawingText())
        parts = [
            '<?xml version="1.0" encoding="UTF-8"?>',
            f'<svg xmlns="http://www.w3.org/2000/svg" width="{total_w}" height="{total_h}">',
            '<rect width="100%" height="100%" fill="white"/>',
            f'<svg x="0" y="0" width="{panel_w}" height="{panel_h}" '
            f'viewBox="{viewbox}" preserveAspectRatio="xMidYMid meet">{content}</svg>',
        ]
        y = panel_h
        if legend_h:
            parts.append(_svg_legend_items(legend_items, LEGEND_MARGIN, y))
            y += legend_h
        if annot_h:
            parts.append(_svg_annotation(annotation, total_w, y))
        parts.append("</svg>")
        with open(path, "w", encoding="utf-8") as f:
            f.write("\n".join(parts))
    else:
        mol_img = Image.open(io.BytesIO(drawer.GetDrawingText())).convert("RGB")
        canvas = Image.new("RGB", (total_w, total_h), "white")
        canvas.paste(mol_img, (0, 0))
        draw = ImageDraw.Draw(canvas)
        y = panel_h
        if legend_h:
            _draw_png_legend(draw, legend_items, total_w, y)
            y += legend_h
        if annot_h:
            _draw_png_annotation(draw, annotation, total_w, y + 4)
        canvas.save(path)

# ---------- reaction helpers ----------
"""Compose a reaction figure from individually-styled molecule panels.

RDKit's built-in reaction drawer (`rdMolDraw2D.DrawReaction`) shares molecule-level
options with `MolDraw2D` (bond length, palette, font), so Post 1's fixes carry over
directly. But it draws catalysts/reagents ("agents") as full structures in the same
row as reactants and products, which crowds the figure fast. This module strips
agents out of the reaction and re-adds them as compact text above the arrow instead,
which is closer to how reactions are drawn in a paper or a slide.

The output format is inferred from the destination filename: ``.svg`` produces a
scalable vector graphic, anything else produces a PNG.
"""
import html
import io
import os
import re
from PIL import Image, ImageDraw, ImageFont
from rdkit import Chem
from rdkit.Chem import AllChem, rdFMCS, rdMolDescriptors

from rdkit.Chem.Draw import rdMolDraw2D

PANEL_SIZE = (260, 220)
ARROW_WIDTH = 140
ARROW_TEXT_PAD = 40
GAP = 16
MARGIN = 20
PLUS_W = 40


def _infer_fmt(path: str, fmt: str | None = None) -> str:
    """Return 'svg' for ``.svg`` paths, otherwise 'png'."""
    if fmt:
        return fmt.lower()
    return "svg" if path.lower().endswith(".svg") else "png"


def _parse_markup(text: str) -> list[tuple[str, str]]:
    """Parse lightweight markup into text segments.

    ``~text~`` becomes a subscript, ``^text^`` becomes a superscript.
    Escaped delimiters ``\\~`` and ``\\^`` are rendered literally.
    Returns a list of ``(segment, style)`` where style is one of
    ``normal``, ``sub``, or ``super``.
    """
    if not text:
        return [("", "normal")]

    segments = []
    pos = 0
    # Match a delimiter that is not escaped by a backslash, then the shortest
    # content up to the next non-escaped matching delimiter.
    pattern = re.compile(r'(?<!\\)([~^])(?P<content>.*?)(?<!\\)\1')
    for match in pattern.finditer(text):
        start, end = match.span()
        if start > pos:
            segments.append((text[pos:start].replace(r'\\', '\\').replace(r'\~', '~').replace(r'\^', '^'), "normal"))
        delim = match.group(1)
        content = match.group("content").replace(r'\\', '\\').replace(r'\~', '~').replace(r'\^', '^')
        segments.append((content, "sub" if delim == "~" else "super"))
        pos = end
    if pos < len(text):
        segments.append((text[pos:].replace(r'\\', '\\').replace(r'\~', '~').replace(r'\^', '^'), "normal"))
    return segments


def _load_font(size):
    try:
        return ImageFont.truetype("DejaVuSans.ttf", size)
    except OSError:
        return ImageFont.load_default()


def _draw_mol_highlighted_bytes(
    mol: Chem.Mol,
    highlight_config: dict,
    fmt: str,
    style: dict | None = None,
) -> bytes:
    """Render an already-abbreviated mol with highlights, returning raw bytes.

    Accepts keys ``atom_colors``, ``bond_colors``, ``atom_maps``, and
    ``annotation`` from ``highlight_config``. Works directly on the mol object so
    no SMILES round-trip occurs (which would expose internal ``*`` abbreviation
    atoms).
    """
    atom_colors = highlight_config.get("atom_colors", {})
    bond_colors = highlight_config.get("bond_colors", {})
    atom_maps = highlight_config.get("atom_maps", {})
    annotation = highlight_config.get("annotation", None)

    highlight_atom_colors = _to_rdkit_color_map(atom_colors, mol)
    highlight_bond_colors = _to_rdkit_color_map(bond_colors, mol)

    highlight_atoms = list(atom_colors.keys()) if atom_colors else []
    bond_indices = (
        _bond_indices_from_atom_pairs(mol, list(bond_colors.keys()))
        if bond_colors else []
    )

    annot_h = _annotation_height(annotation)
    panel_w, panel_h = PANEL_SIZE
    total_h = panel_h + annot_h

    drawer = _make_drawer(PANEL_SIZE, fmt)
    opts = drawer.drawOptions()
    for key, value in {**DEFAULT_STYLE, **(style or {})}.items():
        setattr(opts, key, value)
    opts.useBWAtomPalette()

    if atom_maps:
        _set_atom_map_notes(mol, atom_maps)

    rdMolDraw2D.PrepareAndDrawMolecule(
        drawer,
        mol,
        highlightAtoms=highlight_atoms,
        highlightBonds=bond_indices,
        highlightAtomColors=highlight_atom_colors,
        highlightBondColors=highlight_bond_colors,
    )
    drawer.FinishDrawing()

    raw = drawer.GetDrawingText()
    if fmt == "svg":
        # Return raw RDKit SVG bytes — annotation is handled by the caller
        # (_panel_svg returns it as a separate string so _draw_svg can place it).
        return raw if isinstance(raw, bytes) else raw.encode("utf-8")
    # PNG route
    mol_img = Image.open(io.BytesIO(raw)).convert("RGB")
    if not annotation:
        return raw
    from PIL import ImageDraw as _ImageDraw
    canvas = Image.new("RGB", (panel_w, total_h), "white")
    canvas.paste(mol_img, (0, 0))
    draw = _ImageDraw.Draw(canvas)
    _draw_png_annotation(draw, annotation, panel_w, panel_h + 4)
    buf = io.BytesIO()
    canvas.save(buf, format="PNG")
    return buf.getvalue()


def _apply_panel_abbreviations(mol: Chem.Mol, abbreviations) -> Chem.Mol:
    """Apply abbreviations to a reaction panel molecule.

    Uses ``max_coverage=1.0`` so that fragments such as HOPh are not silently
    skipped on small reactant molecules where the ring would otherwise exceed
    the default 40% coverage limit.
    """
    return apply_abbreviations(mol, abbreviations, max_coverage=1.0)


def _panel_image(
    mol,
    highlight_config: dict | None = None,
    abbreviations=None,
    style: dict | None = None,
) -> Image.Image:
    """Render one molecule as a PIL Image (PNG route)."""
    if abbreviations is not None:
        mol = _apply_panel_abbreviations(mol, abbreviations)
    if highlight_config:
        png_bytes = _draw_mol_highlighted_bytes(
            mol, highlight_config, fmt="png", style=style
        )
        return Image.open(io.BytesIO(png_bytes)).convert("RGB")
    png_bytes = draw_mol_bytes(mol, size=PANEL_SIZE, fmt="png", style=style)
    return Image.open(io.BytesIO(png_bytes)).convert("RGB")


def _panel_svg(
    mol,
    highlight_config: dict | None = None,
    abbreviations=None,
    style: dict | None = None,
) -> tuple[str, str | None]:
    """Render one molecule as a (svg_text, annotation) tuple (SVG route).

    Returns the raw RDKit SVG string and the optional annotation text separately
    so the reaction canvas can position the annotation correctly.
    """
    annotation = (highlight_config or {}).get("annotation", None)
    if abbreviations is not None:
        mol = _apply_panel_abbreviations(mol, abbreviations)
    if highlight_config:
        svg_bytes = _draw_mol_highlighted_bytes(
            mol, highlight_config, fmt="svg", style=style
        )
        return svg_bytes.decode("utf-8"), annotation
    return draw_mol_bytes(mol, size=PANEL_SIZE, fmt="svg", style=style).decode("utf-8"), annotation


def _is_multifragment(mol) -> bool:
    """Return True if the molecule SMILES contains a dot (multiple fragments)."""
    return "." in Chem.MolToSmiles(mol)


def _align_product_to_reference(product_mol, ref_mol, min_mcs_atoms: int = 4) -> bool:
    """Rigid-body align product_mol onto ref_mol using their pairwise MCS.

    Returns True if alignment succeeded, False if it was skipped or failed.
    """
    if _is_multifragment(ref_mol) or _is_multifragment(product_mol):
        return False

    # Both molecules need 2D coordinates before AlignMol can fit them.
    AllChem.Compute2DCoords(ref_mol)
    AllChem.Compute2DCoords(product_mol)

    mcs = rdFMCS.FindMCS(
        [ref_mol, product_mol],
        atomCompare=rdFMCS.AtomCompare.CompareElements,
        bondCompare=rdFMCS.BondCompare.CompareOrderExact,
    )

    if mcs.numAtoms < min_mcs_atoms:
        return False

    mcs_mol = Chem.MolFromSmarts(mcs.smartsString)
    ref_match = ref_mol.GetSubstructMatch(mcs_mol)
    prod_match = product_mol.GetSubstructMatch(mcs_mol)

    if len(ref_match) < mcs.numAtoms or len(prod_match) < mcs.numAtoms:
        return False

    atom_map = list(zip(prod_match, ref_match))
    try:
        AllChem.AlignMol(product_mol, ref_mol, atomMap=atom_map)
        return True
    except Exception:
        return False


def _agent_label(rxn) -> str:
    """Turn agent molecules into a short text label, e.g. 'Pd'."""
    labels = []
    for agent in rxn.GetAgents():
        try:
            Chem.SanitizeMol(agent)
        except Exception:
            pass
        formula = rdMolDescriptors.CalcMolFormula(agent)
        labels.append(formula)
    return ", ".join(labels)


def _layout(n_reactants: int, n_products: int, arrow_width: int = ARROW_WIDTH):
    """Compute the geometry of the final figure.

    Returns ``(total_width, total_height, reactant_xs, product_xs, arrow_x0, arrow_x1, arrow_y)``.
    """
    panel_w, panel_h = PANEL_SIZE
    total_width = (
        MARGIN
        + n_reactants * panel_w + (n_reactants - 1) * PLUS_W
        + GAP + arrow_width + GAP
        + n_products * panel_w + (n_products - 1) * PLUS_W
        + MARGIN
    )
    total_height = panel_h + 2 * MARGIN

    y = MARGIN
    arrow_y = y + panel_h // 2

    reactant_xs = []
    x = MARGIN
    for i in range(n_reactants):
        reactant_xs.append(x)
        x += panel_w
        if i < n_reactants - 1:
            x += PLUS_W

    arrow_x0 = x + GAP
    arrow_x1 = arrow_x0 + arrow_width

    product_xs = []
    x = arrow_x1 + GAP + 10
    for i in range(n_products):
        product_xs.append(x)
        x += panel_w
        if i < n_products - 1:
            x += PLUS_W

    return total_width, total_height, reactant_xs, product_xs, arrow_x0, arrow_x1, arrow_y


def _required_arrow_width(agent_text: str, condition_text: str) -> int:
    """Return an arrow width that fits the widest label plus padding."""
    tmp = Image.new("RGB", (1, 1))
    tmp_draw = ImageDraw.Draw(tmp)
    label_font = _load_font(16)
    max_text_w = 0.0
    if agent_text:
        max_text_w = max(max_text_w, _measure_markup_text(tmp_draw, agent_text, label_font))
    if condition_text:
        for line in condition_text.split("\n"):
            max_text_w = max(max_text_w, _measure_markup_text(tmp_draw, line, label_font))
    return max(ARROW_WIDTH, int(max_text_w + 2 * ARROW_TEXT_PAD))


def _measure_markup_text(
    draw: ImageDraw.ImageDraw,
    text: str,
    base_font: ImageFont.FreeTypeFont,
) -> float:
    """Return total pixel width of a string that may contain ~sub~ or ^super^ markup."""
    sub_super_font = _load_font(int(base_font.size * 0.7))
    total = 0.0
    for segment, style in _parse_markup(text):
        font = sub_super_font if style in ("sub", "super") else base_font
        total += draw.textlength(segment, font=font)
    return total


def _draw_markup_text(
    draw: ImageDraw.ImageDraw,
    text: str,
    x: float,
    y: float,
    base_font: ImageFont.FreeTypeFont,
    fill: str = "black",
) -> None:
    """Draw text with ``~sub~`` and ``^super^`` markup on a PIL canvas."""
    sub_super_font = _load_font(int(base_font.size * 0.7))
    offset = base_font.size * 0.3
    for segment, style in _parse_markup(text):
        font = sub_super_font if style in ("sub", "super") else base_font
        seg_y = y
        if style == "sub":
            seg_y = y + offset
        elif style == "super":
            seg_y = y - offset
        draw.text((x, seg_y), segment, font=font, fill=fill)
        x += draw.textlength(segment, font=font)


def _draw_png(
    path: str,
    reactant_panels: list[Image.Image],
    product_panels: list[Image.Image],
    agent_text: str,
    condition_text: str,
) -> None:
    """Assemble the reaction as a PNG using PIL."""
    n_reactants, n_products = len(reactant_panels), len(product_panels)
    arrow_width = _required_arrow_width(agent_text, condition_text)
    total_width, total_height, reactant_xs, product_xs, arrow_x0, arrow_x1, arrow_y = _layout(
        n_reactants, n_products, arrow_width
    )

    canvas = Image.new("RGB", (total_width, total_height), "white")
    draw = ImageDraw.Draw(canvas)
    label_font = _load_font(16)
    panel_h = PANEL_SIZE[1]
    y = MARGIN

    for x, panel in zip(reactant_xs, reactant_panels):
        canvas.paste(panel, (x, y))

    draw.line([(arrow_x0, arrow_y), (arrow_x1, arrow_y)], fill="black", width=2)
    draw.polygon(
        [(arrow_x1, arrow_y - 6), (arrow_x1, arrow_y + 6), (arrow_x1 + 10, arrow_y)],
        fill="black",
    )

    if agent_text:
        tw = _measure_markup_text(draw, agent_text, label_font)
        x = arrow_x0 + arrow_width / 2 - tw / 2
        _draw_markup_text(draw, agent_text, x, arrow_y - 22, label_font)
    if condition_text:
        for i, line in enumerate(condition_text.split("\n")):
            tw = _measure_markup_text(draw, line, label_font)
            x = arrow_x0 + arrow_width / 2 - tw / 2
            _draw_markup_text(draw, line, x, arrow_y + 22 + i * 18, label_font)

    for x, panel in zip(product_xs, product_panels):
        canvas.paste(panel, (x, y))

    # Plus signs sit halfway between adjacent panels.
    plus_half = 5
    plus_width = 2
    for xs in (reactant_xs, product_xs):
        for i in range(len(xs) - 1):
            cx = xs[i] + PANEL_SIZE[0] + PLUS_W // 2
            cy = y + panel_h // 2
            draw.line(
                [(cx - plus_half, cy), (cx + plus_half, cy)],
                fill="black",
                width=plus_width,
            )
            draw.line(
                [(cx, cy - plus_half), (cx, cy + plus_half)],
                fill="black",
                width=plus_width,
            )

    canvas.save(path)


def _parse_rdkit_svg(svg_text: str) -> tuple[str, str]:
    """Extract the viewBox (or width/height) and inner content from an RDKit molecule SVG."""
    # RDKit may use single or double quotes for SVG attributes.
    viewbox_match = re.search(r"viewBox=(?:'|\")([^'\"]+)(?:'|\")", svg_text)
    if viewbox_match:
        viewbox = viewbox_match.group(1)
    else:
        # Fall back to width/height so the content still scales to fit the panel.
        width_match = re.search(r"width=(?:'|\")(\d+)(?:'|\")", svg_text)
        height_match = re.search(r"height=(?:'|\")(\d+)(?:'|\")", svg_text)
        w = int(width_match.group(1)) if width_match else 100
        h = int(height_match.group(1)) if height_match else 100
        viewbox = f"0 0 {w} {h}"

    # Find the end of the opening <svg ...> tag, then the closing </svg>.
    tag_start = svg_text.find("<svg")
    start = svg_text.find(">", tag_start)
    end = svg_text.rfind("</svg>")
    content = svg_text[start + 1 : end].strip() if start != -1 and end != -1 else ""
    return viewbox, content


def _svg_markup_text(text: str, base_size: int = 16) -> str:
    """Convert ``~sub~`` and ``^super^`` markup into SVG ``<tspan>`` elements."""
    sub_size = int(base_size * 0.7)
    fragments = []
    for segment, style in _parse_markup(text):
        escaped = html.escape(segment)
        if style == "sub":
            fragments.append(f'<tspan baseline-shift="sub" font-size="{sub_size}">{escaped}</tspan>')
        elif style == "super":
            fragments.append(f'<tspan baseline-shift="super" font-size="{sub_size}">{escaped}</tspan>')
        else:
            fragments.append(escaped)
    return "".join(fragments)


def _draw_svg(
    path: str,
    reactant_svgs: list[str],
    product_svgs: list[str],
    agent_text: str,
    condition_text: str,
) -> None:
    """Assemble the reaction as an SVG by inlining each molecule SVG."""
    n_reactants, n_products = len(reactant_svgs), len(product_svgs)
    arrow_width = _required_arrow_width(agent_text, condition_text)
    total_width, total_height, reactant_xs, product_xs, arrow_x0, arrow_x1, arrow_y = _layout(
        n_reactants, n_products, arrow_width
    )

    panel_w, panel_h = PANEL_SIZE
    y = MARGIN

    parts = [
        '<?xml version="1.0" encoding="UTF-8"?>',
        f'<svg xmlns="http://www.w3.org/2000/svg" width="{total_width}" height="{total_height}">',
        '<rect width="100%" height="100%" fill="white"/>',
    ]

    def _embed_panel(x: int, y: int, panel: tuple[str, str | None]):
        svg_text, annotation = panel
        viewbox, content = _parse_rdkit_svg(svg_text)
        result = (
            f'<svg x="{x}" y="{y}" width="{panel_w}" height="{panel_h}" '
            f'viewBox="{viewbox}" preserveAspectRatio="xMidYMid meet">{content}</svg>'
        )
        if annotation:
            annot_y = y + panel_h + ANNOTATION_HEIGHT - 4
            result += (
                f'\n<text x="{x + panel_w // 2}" y="{annot_y}" text-anchor="middle" '
                f'font-size="14" font-family="sans-serif">{html.escape(annotation)}</text>'
            )
        return result

    for x, panel in zip(reactant_xs, reactant_svgs):
        parts.append(_embed_panel(x, y, panel))

    # Arrow shaft and head.
    parts.append(
        f'<line x1="{arrow_x0}" y1="{arrow_y}" x2="{arrow_x1}" y2="{arrow_y}" stroke="black" stroke-width="2"/>'
    )
    parts.append(
        f'<polygon points="{arrow_x1},{arrow_y - 6} {arrow_x1},{arrow_y + 6} {arrow_x1 + 10},{arrow_y}" fill="black"/>'
    )

    if agent_text:
        parts.append(
            f'<text x="{arrow_x0 + arrow_width / 2}" y="{arrow_y - 22}" '
            'text-anchor="middle" font-size="16" font-family="sans-serif">'
            f'{_svg_markup_text(agent_text)}</text>'
        )
    if condition_text:
        for i, line in enumerate(condition_text.split("\n")):
            parts.append(
                f'<text x="{arrow_x0 + arrow_width / 2}" y="{arrow_y + 22 + i * 18}" '
                'text-anchor="middle" font-size="16" font-family="sans-serif">'
                f'{_svg_markup_text(line)}</text>'
            )

    for x, panel in zip(product_xs, product_svgs):
        parts.append(_embed_panel(x, y, panel))

    # Plus signs between panels.
    plus_half = 5
    plus_width = 2
    for xs in (reactant_xs, product_xs):
        for i in range(len(xs) - 1):
            cx = xs[i] + panel_w + PLUS_W // 2
            cy = y + panel_h // 2
            parts.append(
                f'<line x1="{cx - plus_half}" y1="{cy}" x2="{cx + plus_half}" y2="{cy}" '
                f'stroke="black" stroke-width="{plus_width}" stroke-linecap="round"/>'
            )
            parts.append(
                f'<line x1="{cx}" y1="{cy - plus_half}" x2="{cx}" y2="{cy + plus_half}" '
                f'stroke="black" stroke-width="{plus_width}" stroke-linecap="round"/>'
            )

    parts.append("</svg>")

    with open(path, "w", encoding="utf-8") as f:
        f.write("\n".join(parts))


def draw_reaction(
    rxn_smiles: str,
    path: str,
    agent_text: str | None = None,
    condition_text: str | None = None,
    fmt: str | None = None,
    align: bool = True,
    abbreviations=None,
    reactant_abbreviations=None,
    product_abbreviations=None,
    reactant_highlight_configs: list[dict] | None = None,
    product_highlight_configs: list[dict] | None = None,
    style: dict | None = None,
) -> None:
    """Render a reaction as reactant panels -> arrow (with agents/conditions) -> product panels.

    Agents parsed from the reaction SMILES (catalysts, etc.) are shown as a compact
    text label above the arrow rather than as full structures in the main row. If
    ``agent_text`` is supplied, it overrides the auto-derived molecular-formula label.

    The output format is inferred from the file extension unless ``fmt`` is supplied.
    When ``align`` is True, each single-fragment product is rigidly aligned onto the
    first reactant using their maximum common substructure.

    ``abbreviations`` is applied to every panel when supplied.
    ``reactant_abbreviations`` and ``product_abbreviations`` override ``abbreviations``
    for the reactant and product sides respectively, allowing different labels on each
    side (e.g. ``PhOH`` on reactants, ``HOPh`` on the product).
    ``reactant_highlight_configs`` and ``product_highlight_configs`` are lists of
    keyword-argument dicts passed per panel in order. Each config may contain
    ``atom_colors``, ``bond_colors``, ``atom_maps``, and ``annotation``. A
    global ``style`` dict (e.g. ``{"annotationFontScale": 0.7}``) is applied to
    every panel.
    """
    fmt = _infer_fmt(path, fmt)

    rxn = AllChem.ReactionFromSmarts(rxn_smiles, useSmiles=True)
    agent_text = agent_text if agent_text is not None else _agent_label(rxn)

    # Work on copies so alignment never mutates the reaction object itself.
    reactant_mols = [Chem.Mol(m) for m in rxn.GetReactants()]
    product_mols = [Chem.Mol(m) for m in rxn.GetProducts()]

    # Align each product onto the first reactant so the unchanged heavy-atom
    # scaffold stays in the same orientation across the arrow.
    if align and reactant_mols:
        ref_mol = reactant_mols[0]
        for prod_mol in product_mols:
            _align_product_to_reference(prod_mol, ref_mol)

    r_abbrevs = reactant_abbreviations if reactant_abbreviations is not None else abbreviations
    p_abbrevs = product_abbreviations if product_abbreviations is not None else abbreviations
    r_configs = reactant_highlight_configs or [{}] * len(reactant_mols)
    p_configs = product_highlight_configs or [{}] * len(product_mols)

    if fmt == "svg":
        reactant_panels = [
            _panel_svg(m, cfg or None, r_abbrevs, style=style)
            for m, cfg in zip(reactant_mols, r_configs)
        ]
        product_panels = [
            _panel_svg(m, cfg or None, p_abbrevs, style=style)
            for m, cfg in zip(product_mols, p_configs)
        ]
        _draw_svg(path, reactant_panels, product_panels, agent_text, condition_text)
    else:
        reactant_panels = [
            _panel_image(m, cfg or None, r_abbrevs, style=style)
            for m, cfg in zip(reactant_mols, r_configs)
        ]
        product_panels = [
            _panel_image(m, cfg or None, p_abbrevs, style=style)
            for m, cfg in zip(product_mols, p_configs)
        ]
        _draw_png(path, reactant_panels, product_panels, agent_text, condition_text)

Next steps

Part 4 packages the molecule and reaction drawing helpers into a single reusable toolkit. That toolkit will then be released as a free, open-source Python library called acacia_chem on GitHub. Check back soon for the repository and pip install acacia-chem instructions.