autonomous-ai / rdkit

rdkit

RDKit is the cheminformatics toolkit: a molecule is a graph (Chem.Mol), written as SMILES, queried
with SMARTS, given coordinates by distance geometry. Tools: $RDKIT_PYTHON (the pinned venv, with
numpy and pandas), harness_rdkit on PYTHONPATH (build, embed, write, compare),
$RDKIT_TOOLCHAIN/verdict.py (the pane header). Never install another RDKit.

Build, write, verdict

"$RDKIT_PYTHON" molecules/hello.py                        # runs the script → out/<name>.sdf and the rest
"$RDKIT_PYTHON" "$RDKIT_TOOLCHAIN/verdict.py"             # judges the newest molecule → pane header
"$RDKIT_PYTHON" "$RDKIT_TOOLCHAIN/harness_rdkit.py" design "CCO" ethanol   # the same, without a script

from harness_rdkit import design, mol_from_smiles, embed_conformers, embed_3d, write_outputs, properties, similarity, substructure

design("CC(C)Cc1ccc(cc1)C(C)C(=O)O", "ibuprofen")         # parse → 10 conformers → minimise → describe → write
design("CC(C)(O)Cc1ccc(C(C)C(=O)O)cc1", "ibuprofen_oh", parent="ibuprofen")   # an analogue: name its parent

mol = mol_from_smiles("CN1C=NC2=C1C(=O)N(C)C(=O)N2C", "caffeine")   # the steps, when work happens between
confs = embed_conformers(mol, n=10, seed=7)               # ETKDGv3 + MMFF94, deduplicated, lowest first, aligned
write_outputs(confs, "caffeine", parent=False)            # False: not an analogue of anything here

What design/write_outputs write to out/, and what each is for:

The pane shows the molecule while it is being made (out/.progress.json: embedding, minimising,
describing); never write these files by hand. properties(mol) returns the properties without
writing; depict(mol, path) the PNG alone; describe(mol, name) the record and SVG without writing.
read_smi("molecules/series.smi") reads a SMILES name list, table(mols) turns molecules into rows
for pandas.DataFrame.

Parents. An analogue is compared with its parent: the maximum common substructure, the atoms that
changed (+O, −C +O), the property deltas, a 2D depiction drawn on the parent's core and a 3D pose
superposed on it. Pass parent="<name already designed>" (or a SMILES) whenever you make an analogue;
left out, the parent is inferred as the earlier molecule this one is the smallest edit of, and
parent=False says there is none. Design the lead first, then its analogues, so the series reads in
order.

A shape the user kept

The pane can scan a non-ring single bond with native MMFF94, with the other internal coordinates
fixed. The user chooses four connected atoms, moves along the energy curve, and keeps a pose and
note in out/torsions/<id>/. This is a rigid, vacuum, single-bond experiment, not a geometry
optimization or a free-energy calculation. Do not infer populations, kinetics, activity or binding.
The input must have explicit hydrogen atoms, 3D coordinates, 4–200 atoms, one connected component
and MMFF94 parameters. It uses 10°, 15° or 30° spacing and never silently falls back to UFF.

When continuing the user's experiment, read study.json for the method, selected angle, atom
indices, full-precision coordinates and note. Load selected.sdf with
Chem.SDMolSupplier(path, removeHs=False)[0], then author the requested follow-up under a new name.
Keep the original study unchanged. source.mol is the exact input conformer; scan.sdf contains
every sampled angle; energies.csv gives absolute and within-scan relative energies in kcal/mol.
study.zip includes these files and standalone reproduce.py/harness_torsion.py; the recorded
RDKit version is required to verify the calculation. JSON coordinates preserve full precision;
MOL/SDF exchange files round them to four decimals. Unsaved pane scans are not files: ask the user
to Keep study only when their chosen pose or note is needed for the next step.

SMILES, the parts that matter

• Atoms are bare symbols; lowercase is aromatic (c1ccccc1 benzene). Bonds: - single (implicit),
  = double, # triple, / \ around a double bond for E/Z.

• Branches in parentheses: CC(C)C isobutane. Rings close on matching digits: C1CCCCC1 cyclohexane,
  c1ccc2ccccc2c1 naphthalene. Reuse a digit once it is closed.

• Brackets for anything not a default: charge [NH4+], [O-]; isotope [13C]; explicit H [nH] —
  pyrrole is c1cc[nH]c1, and forgetting the H is the commonest SMILES error there is.

• Stereo: [C@H] / [C@@H] at a centre, F/C=C/F trans. Write it when it matters; RDKit will not
  guess, and an unspecified centre silently becomes a racemate.

• Dot separates components: a salt is CC(=O)[O-].[Na+], and most calculations want the parent only.

Useful anchors: water O, ethanol CCO, benzene c1ccccc1, phenol Oc1ccccc1, aspirin
CC(=O)Oc1ccccc1C(=O)O, paracetamol CC(=O)Nc1ccc(O)cc1, caffeine CN1C=NC2=C1C(=O)N(C)C(=O)N2C,
ibuprofen CC(C)Cc1ccc(cc1)C(C)C(=O)O, naproxen COc1ccc2cc(ccc2c1)C(C)C(=O)O, glucose
OC[C@H]1OC(O)[C@H](O)[C@@H](O)[C@@H]1O, penicillin G core CC1(C)S[C@@H]2[C@H](NC(=O)Cc3ccccc3)C(=O)N2[C@H]1C(=O)O.

SMARTS, for finding things

SMARTS is SMILES plus queries: [#6] any carbon, [C,N] either, [!c] not aromatic carbon, [R2]
in two rings, [X3] three connections, [OX2H] a hydroxyl oxygen, * anything, ~ any bond.

substructure(mol, "[OX2H]")                  # hydroxyls → ((3,), (7,))
substructure(mol, "c1ccccc1")                # benzene rings
substructure(mol, "[CX3](=O)[OX2H1]")        # carboxylic acid

Groups worth keeping: carboxylic acid [CX3](=O)[OX2H1], amide [NX3][CX3](=[OX1]), primary amine
[NX3;H2;!$(NC=O)], sulfonamide [SX4](=[OX1])(=[OX1])([NX3]), nitro [N+](=O)[O-], halogen [F,Cl,Br,I].

Common tasks

Modify a scaffold — edit the SMILES where the substituent goes, or replace a group in place:

from rdkit import Chem
core = Chem.MolFromSmiles("CC(=O)Oc1ccccc1C(=O)O")
out  = Chem.ReplaceSubstructs(core, Chem.MolFromSmarts("[CX3](=O)[OX2H1]"),
                              Chem.MolFromSmiles("C(=O)NC"), replaceAll=True)[0]
Chem.SanitizeMol(out); print(Chem.MolToSmiles(out))

Enumerate analogues — one substituent list, one loop, a table, and the ones worth looking at written
out, each with its parent:

import pandas as pd
from rdkit import Chem
from harness_rdkit import mol_from_smiles, properties, design
subs = {"H": "", "F": "F", "Cl": "Cl", "OMe": "OC", "CF3": "C(F)(F)F"}
mols = [mol_from_smiles(f"CC(C)Cc1ccc(C(C)C(=O)O)c({r})c1" if r else "CC(C)Cc1ccc(cc1)C(C)C(=O)O", name)
        for name, r in subs.items()]
print(pd.DataFrame([{"R": m.GetProp("_Name"), **properties(m)} for m in mols])[["R", "mw", "logp", "tpsa", "qed"]])
design("CC(C)Cc1ccc(cc1)C(C)C(=O)O", "ibuprofen")         # the lead first
design(Chem.MolToSmiles(mols[-1]), "analogue_cf3", parent="ibuprofen")   # then the analogue, into the series

Similarity search across a list — Morgan (ECFP4) Tanimoto; > 0.7 is a close analogue, < 0.3 unrelated:

from harness_rdkit import read_smi, similarity
query = "CC(=O)Oc1ccccc1C(=O)O"
hits = sorted(((similarity(query, m), m.GetProp("_Name")) for m in read_smi("molecules/library.smi")), reverse=True)
for score, name in hits[:10]: print(f"{score:.2f}  {name}")

3D and conformers — design already runs a small conformer search (conformers=10): ETKDGv3
embeddings, MMFF94 minimisation, duplicates dropped, the lowest within 10 kcal/mol kept, aligned, and
written as the ensemble the pane plays. Raise it for a flexible molecule (conformers=30), lower it to
1 for a single seed-7 pose. embed_conformers(mol, n, seed) is the same search by hand; an
unspecified stereocentre is pinned to one configuration across the ensemble and still reported as
unspecified. Energies are vacuum force-field numbers — a guide to shape, not to populations.

Save for other tools: write_outputs writes the SDFs. Chem.MolToPDBFile(mol, "out/x.pdb"),
Chem.MolToXYZFile(mol, "out/x.xyz"), Chem.MolToSmiles(mol) for the canonical string.

Pitfalls

• Sanitization is not optional. Chem.MolFromSmiles returns None for anything that will not
  sanitize — a five-valent carbon, an unclosed ring, c1ccccc1 written with the wrong aromaticity.
  A None is a typo in the SMILES; fix the string, never sanitize=False your way past it.

• Hydrogens. Properties are computed on the graph without them; 3D needs them. embed_3d adds
  them and gives you a new molecule — the original stays flat, which is what the depiction wants.

• Stereo. An unspecified centre embeds as one arbitrary enantiomer and the SDF will look decided.
  Write the stereo into the SMILES, or say in your answer that the centre is unspecified.

• Protonation. SMILES are written neutral by convention; a carboxylic acid is C(=O)O even though
  it is an anion at pH 7.4. cLogP and TPSA assume the neutral form — say so rather than "correcting" it.

• Salts and mixtures. Strip them before computing: Chem.MolStandardize.rdMolStandardize.LargestFragmentChooser().

• Macrocycles and cages sometimes fail to embed; embed_3d retries with random coordinates, and
  after that it is a different seed or useMacrocycleTorsions.

• cLogP is Crippen's estimate, QED a 2012 desirability score, Lipinski a rule of thumb for oral
  absorption. They are guides, not measurements, and nothing here predicts activity, binding or safety.

Rules

• molecules/ holds scripts and .smi lists, out/ holds everything produced. One design call per
  molecule, name matching the file you want in the pane.

• Run the verdict after every script: "$RDKIT_PYTHON" "$RDKIT_TOOLCHAIN/verdict.py". It is ready when
  a script exists, the newest out/*.sdf parses with a real 3D conformer, and no error is open;
  Lipinski and Veber violations, PAINS/Brenk alerts and a strained energy are warnings, an unspecified
  stereocentre a note — none of them failures.

• The pane opens the newest molecule and follows each new one: 3D (wire, stick, ball-and-stick,
  spacefill; hydrogens; colour by element, Gasteiger charge, Crippen lipophilicity or change from the
  parent; VDW/SAS/SES surfaces coloured by electrostatic or lipophilicity potential; hover for an atom's
  charge and hybridisation; click to measure distances, angles and dihedrals), the 2D depiction linked
  atom-for-atom with its functional groups, the properties with Lipinski/Veber badges and plain-language
  flags, the conformer player, the parent overlay, and the series as a strip and a comparison table.
  Writing the files through design/write_outputs is how you put something in it; point the user at
  what to look at ("colour by change", "the conformer tab") rather than describing the picture.