Enumerating Atom Substitutions

Problem

You want to enumerate all possible carbon-to-nitrogen atom substitutions in a molecule. See examples in Table 1.

Table 1. Example of enumerating atom substitutions

one substitution

two substitutions

three substitutions

four substitutions
../_images/enumsubstitutions-one.svg ../_images/enumsubstitutions-two.svg ../_images/enumsubstitutions-three.svg ../_images/enumsubstitutions-four.svg

Ingredients

Difficulty Level

🌶️

Download

Download code

enumsubstitutions.py

See also Usage subsection.

Source Code

enumsubstitutions
#!/usr/bin/env python3
# (C) 2026 Cadence Design Systems, Inc. (Cadence)
# All rights reserved.
# TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
# provided to current licensees or subscribers of Cadence products or
# SaaS offerings (each a "Customer").
# Customer is hereby permitted to use, copy, and modify the Sample Code,
# subject to these terms. Cadence claims no rights to Customer's
# modifications. Modification of Sample Code is at Customer's sole and
# exclusive risk. Sample Code may require Customer to have a then
# current license or subscription to the applicable Cadence offering.
# THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
# NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
# PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
# liable for any damages or liability in connection with the Sample Code
# or its use.

"""Enumerate carbon-to-nitrogen substitutions in a molecule."""

import argparse
import os
import pathlib
import sys
from itertools import combinations

from openeye import oechem
from rich.console import Console
from rich_argparse import HelpPreviewAction, RichHelpFormatter

__SCRIPT_NAME__ = pathlib.Path(__file__).absolute().stem
__SCRIPT_DESC__ = "Enumerate carbon-to-nitrogen substitutions in a molecule."
__SCRIPT_TOOLKITS__ = ["oechem"]
__SCRIPT_CATEGORIES__ = ["cheminformatics"]


def parse_options() -> argparse.Namespace:
    """Set up command line options."""
    parser = argparse.ArgumentParser(
        add_help=True,
        formatter_class=RichHelpFormatter,
        description="[yellow]" + __SCRIPT_DESC__ + "[/yellow]",
    )

    io_group = parser.add_argument_group("Input/output options")
    io_group.add_argument(
        "--in",
        dest="input",
        metavar="INPUT-FILE",
        type=str,
        required=True,
        help="input molecule file",
    )
    io_group.add_argument(
        "--out",
        dest="output",
        metavar="OUTPUT-FILE",
        type=str,
        required=True,
        help="output molecule file",
    )

    sub_group = parser.add_argument_group("Substitution options")
    sub_group.add_argument(
        "--min-subs",
        metavar="N",
        type=int,
        required=True,
        help="minimum number of carbon-to-nitrogen substitutions",
    )
    sub_group.add_argument(
        "--max-subs",
        metavar="N",
        type=int,
        required=True,
        help="maximum number of carbon-to-nitrogen substitutions",
    )

    parser.add_argument("--help-image", action=HelpPreviewAction)
    parser.add_argument(
        "--save-console-svg",
        default=False,
        action="store_true",
        help=f"run command and capture console output in {__SCRIPT_NAME__}.svg file",
    )
    return parser.parse_args()


def enumerate_substitutions(
    mol: oechem.OEMolBase,
    min_subs: int,
    max_subs: int,
    atom_pred: oechem.OEUnaryAtomPred | None = None,
) -> list[oechem.OEMolBase]:
    """Enumerate all unique carbon-to-nitrogen substitution combinations."""
    if atom_pred is None:
        atom_pred = oechem.OEIsCarbon()
    oechem.OESuppressHydrogens(mol)
    mol.Sweep()
    atom_indices = [atom.GetIdx() for atom in mol.GetAtoms(atom_pred)]

    substituted_mols: list[oechem.OEMolBase] = []
    for n in range(min_subs, max_subs + 1):
        for atom_comb in combinations(atom_indices, n):
            sub_mol = oechem.OEGraphMol(mol)
            if perform_substitution(sub_mol, atom_comb):
                substituted_mols.append(oechem.OEGraphMol(sub_mol))

    return substituted_mols


def perform_substitution(mol: oechem.OEMolBase, atom_comb: tuple[int, ...]) -> bool:
    """Substitute selected carbon atoms with nitrogen."""
    nitrogen_valence = 3

    for atom_idx in atom_comb:
        atom = mol.GetAtom(oechem.OEHasAtomIdx(atom_idx))
        if atom is None:
            return False
        atom.SetAtomicNum(oechem.OEElemNo_N)

        explicit_val = atom.GetExplicitValence()
        if nitrogen_valence > explicit_val:
            atom.SetImplicitHCount(nitrogen_valence - explicit_val)
            atom.SetFormalCharge(0)
        elif nitrogen_valence < explicit_val:
            atom.SetImplicitHCount(0)
            atom.SetFormalCharge(explicit_val - nitrogen_valence)
        else:
            atom.SetImplicitHCount(0)
            atom.SetFormalCharge(0)

    return True


def write_unique_molecules(
    mols: list[oechem.OEMolBase], ofs: oechem.oemolostream
) -> int:
    """Write only unique molecules (by canonical SMILES) to the output stream."""
    unique_smiles: set[str] = set()
    count = 0
    for mol in mols:
        smi = oechem.OEMolToSmiles(mol)
        if smi not in unique_smiles:
            unique_smiles.add(smi)
            oechem.OEWriteMolecule(ofs, mol)
            count += 1
    return count


def main() -> int:
    """Enumerate carbon-to-nitrogen substitutions."""
    args = parse_options()

    console = Console(record=args.save_console_svg)

    ifs = oechem.oemolistream()
    if not ifs.open(args.input):
        oechem.OEThrow.Fatal(f"Cannot open input file: {args.input}")

    ofs = oechem.oemolostream()
    if not ofs.open(args.output):
        oechem.OEThrow.Fatal(f"Cannot open output file: {args.output}")

    mol = oechem.OEGraphMol()
    if not oechem.OEReadMolecule(ifs, mol):
        oechem.OEThrow.Fatal(f"Cannot read molecule from: {args.input}")

    min_subs = max(1, min(args.min_subs, args.max_subs))
    max_subs = min(mol.NumAtoms(), max(args.min_subs, args.max_subs))

    substituted_mols = enumerate_substitutions(mol, min_subs, max_subs)
    count = write_unique_molecules(substituted_mols, ofs)

    console.print(f"Generated [green]{count}[/green] unique substituted molecules.")

    if args.save_console_svg:
        console.save_svg(f"{__SCRIPT_NAME__}.svg", title=__SCRIPT_NAME__)

    return os.EX_OK


setattr(main, "__SCRIPT_NAME__", __SCRIPT_NAME__)
setattr(main, "__SCRIPT_DESC__", __SCRIPT_DESC__)
setattr(main, "__SCRIPT_TOOLKITS__", __SCRIPT_TOOLKITS__)
setattr(main, "__SCRIPT_CATEGORIES__", __SCRIPT_CATEGORIES__)

if __name__ == "__main__":
    sys.exit(main())

Solution

The enumerate_substitutions function generates all possible combinations of the given atom indices in the range of [min_subs, max_subs].

def enumerate_substitutions(
    mol: oechem.OEMolBase,
    min_subs: int,
    max_subs: int,
    atom_pred: oechem.OEUnaryAtomPred | None = None,
) -> list[oechem.OEMolBase]:
    """Enumerate all unique carbon-to-nitrogen substitution combinations."""
    if atom_pred is None:
        atom_pred = oechem.OEIsCarbon()
    oechem.OESuppressHydrogens(mol)
    mol.Sweep()
    atom_indices = [atom.GetIdx() for atom in mol.GetAtoms(atom_pred)]

    substituted_mols: list[oechem.OEMolBase] = []
    for n in range(min_subs, max_subs + 1):
        for atom_comb in combinations(atom_indices, n):
            sub_mol = oechem.OEGraphMol(mol)
            if perform_substitution(sub_mol, atom_comb):
                substituted_mols.append(oechem.OEGraphMol(sub_mol))

    return substituted_mols

The perform_substitution function carries out the carbon-to-nitrogen atom substitutions defined by the given atom combination. After setting the atomic number, the implicit hydrogen count and the charge of the altered atoms are corrected to generate chemically correct molecules.

def perform_substitution(mol: oechem.OEMolBase, atom_comb: tuple[int, ...]) -> bool:
    """Substitute selected carbon atoms with nitrogen."""
    nitrogen_valence = 3

    for atom_idx in atom_comb:
        atom = mol.GetAtom(oechem.OEHasAtomIdx(atom_idx))
        if atom is None:
            return False
        atom.SetAtomicNum(oechem.OEElemNo_N)

        explicit_val = atom.GetExplicitValence()
        if nitrogen_valence > explicit_val:
            atom.SetImplicitHCount(nitrogen_valence - explicit_val)
            atom.SetFormalCharge(0)
        elif nitrogen_valence < explicit_val:
            atom.SetImplicitHCount(0)
            atom.SetFormalCharge(explicit_val - nitrogen_valence)
        else:
            atom.SetImplicitHCount(0)
            atom.SetFormalCharge(0)

    return True

After enumerating all possible atom substitutions, the molecules are written to the output stream by calling the write_unique_molecules function. To avoid outputting duplicate molecules, the OEMolToSmiles function is used to generate canonical isomeric SMILES. The OEMolToSmiles function generates the same SMILES string for isomorphic molecules, i.e. molecules that only differ by atom ordering. See example of graph isomorphic molecules in Table 2.

def write_unique_molecules(
    mols: list[oechem.OEMolBase], ofs: oechem.oemolostream
) -> int:
    """Write only unique molecules (by canonical SMILES) to the output stream."""
    unique_smiles: set[str] = set()
    count = 0
    for mol in mols:
        smi = oechem.OEMolToSmiles(mol)
        if smi not in unique_smiles:
            unique_smiles.add(smi)
            oechem.OEWriteMolecule(ofs, mol)
            count += 1
    return count
Table 2. Example of graph isomorphic molecules.
../_images/substitution-isomorph-A.svg ../_images/substitution-isomorph-B.svg

Usage

See Download section to download the script.

> enumsubstitutions --help
../_images/enumsubstitutions-help.svg

The enumsubstitutions script enumerates all possible carbon-to-nitrogen atom substitutions in a molecule within a specified range of substitution counts.

> enumsubstitutions --in molecule.ism --out output.ism --min-subs 1 --max-subs 2

The enumsubstitutions script reports the number of unique substituted molecules generated.

../_images/enumsubstitutions-01-stdout.svg

The command above generates output.ism:

C1CNCC1=O
C1CC(=O)NC1
C1CC[N+](=O)C1
C1CNNC1=O
C1C[N+](=O)CN1
C1C(=O)NCN1
C1C(=O)CNN1
C1CN[N+](=O)C1
C1CNC(=O)N1

Discussion

The enumerate_substitutions function takes an atom predicate that determines which atoms are substituted. By default, it only allows substitution of carbon atoms. OEChem TK has a diverse selection of built-in atom predicates and also allows the creation of user-defined ones. This provides a flexible and convenient way to determine the set of atoms to be substituted. For example, when using the OEIsHeavy predicate, the enumeration of any two substitutions results in the following extra three molecules for the cyclopentanone input structure:

Table 3. Extra substitutions when allowing substitution of any heavy atom.
../_images/enumsubstitutions-extra-1.svg ../_images/enumsubstitutions-extra-2.svg ../_images/enumsubstitutions-extra-3.svg

See also in Python documentation

  • itertools.combinations()

See also in OEChem TK manual

Theory

API