Skip to content

Add determine_bond_orders kwarg to Sire-RDKit conversion function. #465

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
lohedges merged 1 commit into from
Aug 26, 2026
Merged

Add determine_bond_orders kwarg to Sire-RDKit conversion function. #465

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions doc/source/changelog.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,11 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
* Fixed a dangling ``else`` statement that caused ``setAmberWater`` and
``setGromacsWater`` to fail when using the TIP4P water model.

* Added a ``determine_bond_orders`` keyword argument (and equivalent property map
option) to the Sire-to-RDKit conversion functions. This defaults to ``True``,
but can be set to ``False`` to fall back on the internal bond inference
heuristic, which is much faster for large molecules, e.g. proteins.

`2026.1.0 <https://github.com/openbiosim/sire/compare/2025.4.0...2026.1.0>`__ - June 2026
-----------------------------------------------------------------------------------------

Expand Down
50 changes: 42 additions & 8 deletions src/sire/convert/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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":
Expand All @@ -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."
)


Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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)

Expand All @@ -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
Expand Down
123 changes: 70 additions & 53 deletions wrapper/Convert/SireRDKit/sire_rdkit.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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<SireUnits::Dimension::Charge>(map["charge"]).to(SireUnits::mod_electron);
total_charge += a->getFormalCharge();
}
total_charge = static_cast<int>(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<SireUnits::Dimension::Charge>(map["charge"]).to(SireUnits::mod_electron);
}
total_charge = static_cast<int>(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())
{
Expand All @@ -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 (...)
{
}
}
}

Expand Down
Loading