Ligands: from atoms to SMILES, the results

Forecast: FORECAST.md, written on 2026-10-16 at 18:00, before any ligand was read and left unedited. Code: qo_chem/ligands.py. Scripts: tools/ligands/benchmark.py (the test set) and tools/ligands/heldout.py (the held-out set). Tests: tests/test_ligands.py.

#The question

A ligand in a structure file is a three-letter code and heavy atoms with coordinates. Can the app turn that into a molecule with bond orders, hydrogens and stereo? It needs one to open the ligand on the other pages, which start from a SMILES, and to run the orbital layer where the ligand sits.

#The test set, and what counts as right

• Test set: 104 ligand instances of 24 components from 20 structures. The files come from the test data of MDAnalysis, Biopython, ProDy and OpenMM (tests/data/ligands/README.md). The reference is each component's entry in the PDB's Chemical Component Dictionary (CCD), as SMILES.
• Same molecule: the same heavy-atom graph and the same number of hydrogens on every atom once both are neutralised. Bond orders then agree up to resonance.
• Complete instance: an instance can be scored only when it has all its atoms. Atoms missing at a covalent bond to the protein are the dictionary's leaving atoms; they are allowed (see "What changed", step 1).

#Run 1: the rules as written in the forecast

All five failures were the same mistake: an amino group on an aromatic ring read as an imine. This hit adenine in SAH and in the three ATPγS, and the aminotriazine of ZM241385. The bond lengths cannot tell the two tautomers apart, and the matching took the imine.

Two things the forecast did not foresee:

• 27 instances could not be scored, because atoms were missing. In 21 of them the missing atom sits where the ligand is bonded to the protein: O1 of a sugar on an Asn, or the OH of 4-guanidinobenzoate on Ser195 of trypsin (an acyl-enzyme).
• All 18 haems in the test set are haem c. Each one is bonded to two cysteines through its former vinyl carbons (S–C 1.82–1.87 Å), which run 1 could not see. It gave nine of the 18 both vinyls and six of them one.

#What changed, in order

1. Bonds to the rest of the structure. Atoms of other residues within bonding distance, or declared in LINK / _struct_conn, now count as neighbours when the geometry is read. In the SMILES they become H. A dictionary atom missing at such a link is a leaving atom; one missing elsewhere means the model is incomplete, and the report says so. 37 test instances have a link: 18 haems, 17 NAG, 1 BMA and 1 GBS.
2. An exocyclic N on an aromatic ring is first tried as an amine. This fixes adenine and the aminotriazine. It is the tautomer found in water, and a chemical rule, not a fitted one.

Run 2 (test set):

• 78 of 78 complete instances are right.
• 7 instances are named incomplete: two sulfates that are a lone S, and five NAG with O1 absent and no bond modelled to the Asn.
• All 151 CCD stereocentres have the handedness found in the file, and every written SMILES reads back to the same molecule and stereo.

These numbers are in-sample now: rule 2 was made after seeing the failures. So a second set was taken that the rules had not seen.

#The held-out set

• What it is: 156 instances of 31 components from the 23 test structures of biotite (tests/data/ligands/heldout.json). These files carry the PDB's own bond table for each component (_chem_comp_bond), which is the reference.
• Overlap: five components (NAG, ACT, SO4, BMA, GOL) also occur in the test set, as other instances.

Run 1 of the held-out set (rules as after step 2, before any look at this set):

• complete instances: 153/156; the other 3 are fragments (a γ-phosphate missing from GTP, half a PEG, the adenosyl part of cobalamin);
• components new to the rules: 87/92 (94.6 %);
• all: 148/153 (96.7 %).

The five failures, and why:

Then step 3:

3. Three chemical rules:
   • C=O is preferred over C=N or C=C next to an O–H (amide over imidic acid, lactam over lactim, ketone over enol);
   • no double bond is placed where a bond is longer than a single bond;
   • a ring counts as aromatic only if none of its bonds is that long;
   • a squeezed carbon can still take a double bond that is clearly short.

After step 3: held-out 151/153 (98.7 %); components new to the rules 90/92. FMN, whose deposited geometry really is the lactim, and the corrin remain wrong.

#Haem (forecast L4)

The forecast expected haem b: two N⁻ toward the iron, two imine N, C₃₄H₃₂N₄O₄²⁻. Both halves were wrong.

The test set holds no haem b. In 19HC every haem is haem c. Read with its links, it is C₃₄H₃₆N₄O₄²⁻ once each cysteine is replaced by H: mesoporphyrin, with both former vinyls saturated. The deposited coordinates still keep a vinyl bond length (CAB–CBB 1.37 Å) next to the S–C bond. The model was refined with haem b's restraints plus two links. The report says this in a note.

The porphyrin core. The largest-weight matching gave four C=N (a 2e-oxidised porphyrin, 12 double bonds). A penalty on N double bonds in five-membered rings gave none at all (isophlorin, 10). No linear weight lands on two. So a four-N core around one centre (porphyrin, chlorin, phthalocyanine) now gets exactly two C=N. Every pair is tried and the best kept.

Result: 18/18 haems as dianions with two N⁻ and two C=N. A free-base porphyrin, read without its iron, gets two N–H (test_free_base_porphyrin_has_two_nh). Corrin, which carries one negative charge, is not covered: that is the B1Z failure.

#The forecasts, scored

#What this means for the app

• The order of the routes. The file's own dictionary comes first, then the CCD SMILES, then the name. The first that has the same atoms as the file supplies the bond orders and hydrogens: the dictionary knows the tautomer, and the geometry does not.
• Disagreements are shown. A dictionary that differs from the geometry only next to a covalent link is not used there, because the link changed the chemistry: haem c is not haem b.
• When no dictionary has the same atoms (a wrong code, a "LIG" from a modelling program, an incomplete model), the geometry route is used, and the report says why.
• The free molecule. A dictionary's free molecule, with its leaving atoms back, opens on the other pages. With a metal, its organic part, neutral: protoporphyrin IX for HEM.

#Not verified here

• The real responses of data.rcsb.org (/rest/v1/core/chemcomp/{id}) and of GL1X's ligand index (/api/bio/ligands). The browser test used stand-ins in the format the code expects. The page accepts several field names for the index.
• OPSIN on the server (forecast L7).
• The real Mol* and 3Dmol.js bundles.

The first-run checklist is in DEPLOYMENT.md.