From 4cbe35fb076eb3dcddff4d8bb551981ec1f07327 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 27 Aug 2026 19:11:55 +0100 Subject: [PATCH] Backport fix from PR #465. [ci skip] --- src/sire/convert/__init__.py | 50 +++++++-- wrapper/Convert/SireRDKit/sire_rdkit.cpp | 123 +++++++++++++---------- 2 files changed, 112 insertions(+), 61 deletions(-) diff --git a/src/sire/convert/__init__.py b/src/sire/convert/__init__.py index 243c512a6..e036d9d05 100644 --- a/src/sire/convert/__init__.py +++ b/src/sire/convert/__init__.py @@ -48,7 +48,7 @@ def supported_formats(): return _supported_formats() -def to(obj, format: str = "sire", map=None): +def to(obj, format: str = "sire", map=None, determine_bond_orders: bool = True): """ Convert the passed object from its current object format to the specified object format (default "sire"). Typically this will be converting @@ -62,13 +62,19 @@ def to(obj, format: str = "sire", map=None): The format to convert to map: The property map to use for the conversion + determine_bond_orders: bool (default True) + Whether to use RDKit's ``determineBondOrders`` function when bond + orders need to be inferred during conversion to rdkit format. This + is more robust than the internal heuristic, but can be slow for + large molecules, e.g. proteins. (Only used when converting to + rdkit format.) """ format = format.lower() if format == "sire": return to_sire(obj, map=map) elif format == "rdkit": - return to_rdkit(obj, map=map) + return to_rdkit(obj, map=map, determine_bond_orders=determine_bond_orders) elif format == "gemmi": return to_gemmi(obj, map=map) elif format == "biosimspace": @@ -77,7 +83,7 @@ def to(obj, format: str = "sire", map=None): return to_openmm(obj, map=map) else: raise ValueError( - f"Cannot convert {obj} as the format '{format}' is " "not recognised." + f"Cannot convert {obj} as the format '{format}' is not recognised." ) @@ -184,12 +190,24 @@ def to_biosimspace(obj, map=None): return sire_to_biosimspace(to_sire(obj, map=map), map=map) -def to_rdkit(obj, map=None): +def to_rdkit(obj, map=None, determine_bond_orders: bool = True): """ Convert the passed object from its current object format to a rdkit object format. + + Args: + obj: + The input object to convert + map: + The property map to use for the conversion + determine_bond_orders: bool (default True) + Whether to use RDKit's ``determineBondOrders`` function when bond + orders need to be inferred. This is more robust than the internal + heuristic, but can be slow for large molecules, e.g. proteins. """ - return sire_to_rdkit(to_sire(obj, map=map), map=map) + return sire_to_rdkit( + to_sire(obj, map=map), map=map, determine_bond_orders=determine_bond_orders + ) def to_gemmi(obj, map=None): @@ -425,10 +443,20 @@ def rdkit_to_sire(obj, map=None): return mols -def sire_to_rdkit(obj, map=None): +def sire_to_rdkit(obj, map=None, determine_bond_orders: bool = True): """ Convert the passed sire object (either a molecule or list of molecules) to a rdkit equivalent + + Args: + obj: + The sire object to convert + map: + The property map to use for the conversion + determine_bond_orders: bool (default True) + Whether to use RDKit's ``determineBondOrders`` function when bond + orders need to be inferred. This is more robust than the internal + heuristic, but can be slow for large molecules, e.g. proteins. """ obj = _to_selectormol(obj) @@ -441,9 +469,15 @@ def sire_to_rdkit(obj, map=None): "'conda install -c conda-forge rdkit'" ) - from ..base import create_map + from ..base import create_map, wrap + + map = create_map(map) + + # only use the kwarg if this hasn't already been set in the property map + if not map.specified("determine_bond_orders"): + map.set("determine_bond_orders", wrap(determine_bond_orders)) - mols = _sire_to_rdkit(obj, map=create_map(map)) + mols = _sire_to_rdkit(obj, map=map) if mols is None: return None diff --git a/wrapper/Convert/SireRDKit/sire_rdkit.cpp b/wrapper/Convert/SireRDKit/sire_rdkit.cpp index 27ca5ac99..7e081357f 100644 --- a/wrapper/Convert/SireRDKit/sire_rdkit.cpp +++ b/wrapper/Convert/SireRDKit/sire_rdkit.cpp @@ -678,6 +678,15 @@ namespace SireRDKit force_stereo_inference = map["force_stereo_inference"].value().asABoolean(); } + // Whether to use RDKit's determineBondOrders() to infer bond orders. This + // is more robust than our heuristic, but can be prohibitively slow for + // large molecules, e.g. proteins. + bool determine_bond_orders = true; + if (map.specified("determine_bond_orders")) + { + determine_bond_orders = map["determine_bond_orders"].value().asABoolean(); + } + for (int i = 0; i < atoms.count(); ++i) { const auto atom = atoms(i); @@ -860,34 +869,39 @@ namespace SireRDKit // integer formal charge of the molecule). int total_charge = 0; - if (has_bond_info and force_stereo_inference) + if (determine_bond_orders) { - for (auto a : molecule.atoms()) + if (has_bond_info and force_stereo_inference) { - total_charge += a->getFormalCharge(); - } - } - else - { - try - { - double charge_sum = 0.0; - for (int i = 0; i < atoms.count(); ++i) + for (auto a : molecule.atoms()) { - charge_sum += atoms(i).property(map["charge"]).to(SireUnits::mod_electron); + total_charge += a->getFormalCharge(); } - total_charge = static_cast(std::round(charge_sum)); } - catch (...) + else { - total_charge = 0; + try + { + double charge_sum = 0.0; + for (int i = 0; i < atoms.count(); ++i) + { + charge_sum += atoms(i).property(map["charge"]).to(SireUnits::mod_electron); + } + total_charge = static_cast(std::round(charge_sum)); + } + catch (...) + { + total_charge = 0; + } } } // When bond info is present but force_stereo_inference is requested, - // reset all bonds to SINGLE and clear formal charges so that the - // inference algorithm starts from a clean connectivity graph. - if (has_bond_info and force_stereo_inference) + // reset all bonds to SINGLE and clear formal charges so that + // determineBondOrders() starts from a clean connectivity graph. The + // heuristic in infer_bond_info() works from the unpaired electron + // count of each atom, so doesn't need (or want) this reset. + if (determine_bond_orders and has_bond_info and force_stereo_inference) { for (auto b : molecule.bonds()) { @@ -903,49 +917,52 @@ namespace SireRDKit molecule.updatePropertyCache(false); } - // Prefer RDKit's determineBondOrders, which is based on the xyz2mol - // linear-programming algorithm and is significantly more robust than the - // MDAnalysis heuristic implemented in infer_bond_info(). - // - // determineBondOrders() needs all heavy atoms to have noImplicit set so - // that it does not try to add implicit hydrogens (all H are explicit when - // loaded from formats such as AMBER that carry all hydrogen atoms). - for (auto a : molecule.atoms()) - { - if (a->getAtomicNum() > 1) - { - a->setNoImplicit(true); - } - } + bool inferred = false; - // Check for dummy atoms (atomic_num == 0): determineBondOrders may not - // handle them correctly, so fall back to the heuristic in that case. - bool has_dummy_atoms = false; - for (auto a : molecule.atoms()) + if (determine_bond_orders) { - if (a->getAtomicNum() == 0) + // Prefer RDKit's determineBondOrders, which is based on the xyz2mol + // linear-programming algorithm and is significantly more robust than the + // MDAnalysis heuristic implemented in infer_bond_info(). + // + // determineBondOrders() needs all heavy atoms to have noImplicit set so + // that it does not try to add implicit hydrogens (all H are explicit when + // loaded from formats such as AMBER that carry all hydrogen atoms). + for (auto a : molecule.atoms()) { - has_dummy_atoms = true; - break; + if (a->getAtomicNum() > 1) + { + a->setNoImplicit(true); + } } - } - - bool inferred = false; - if (not has_dummy_atoms) - { - try + // Check for dummy atoms (atomic_num == 0): determineBondOrders may not + // handle them correctly, so fall back to the heuristic in that case. + bool has_dummy_atoms = false; + for (auto a : molecule.atoms()) { - // embedChiral=false: we call sanitizeMol ourselves below, - // and assignStereochemistryFrom3D is called afterwards. - RDKit::determineBondOrders(molecule, total_charge, - /*allowChargedFragments=*/true, - /*embedChiral=*/false, - /*useAtomMap=*/false); - inferred = true; + if (a->getAtomicNum() == 0) + { + has_dummy_atoms = true; + break; + } } - catch (...) + + if (not has_dummy_atoms) { + try + { + // embedChiral=false: we call sanitizeMol ourselves below, + // and assignStereochemistryFrom3D is called afterwards. + RDKit::determineBondOrders(molecule, total_charge, + /*allowChargedFragments=*/true, + /*embedChiral=*/false, + /*useAtomMap=*/false); + inferred = true; + } + catch (...) + { + } } }