Use `rdkit` for SSSR and RCs (bug fix + Python upgrade) by JacksonBurns · Pull Request #2796 · ReactionMechanismGenerator/RMG-Py

JacksonBurns · 2025-05-25T21:43:19Z

Currently we use RingDecomposerLib for finding the Smallest Set of Smallest Rings and getting the Relevant Cycles. This package does not support Python 3.10+ and is thus blocking further upgrades to RMG.

@KnathanM in particular is looking to get RMG to Python 3.11 so as to add support for ChemProp v2.

I believe we can just use RDKit to do these operations instead. The original paper mentions that the functionality was being moved upstream to RDKit. With the help of AI I've taken just a first pass at reimplementing, with the special note that:

I opted to use the Symmetric SSSR in place of the 'true' SSSR. This is because the latter is non-unique (see RDKit's "The SSSR Problem"). This should actually resolve get_deterministic_sssr is not really deterministic #2562
I need to read more about the "Relevant Cycles"

This PR will be a draft for now, as it is predicated on Python 3.9 already being available (which it nearly is in #2741)

Motivation or Problem

A clear and concise description of what what you're trying to fix or improve. Please reference any issues that this addresses.

Description of Changes

A clear and concise description of what what you've changed or added.

Testing

A clear and concise description of testing that you've done or plan to do.

Reviewer Tips

Suggestions for verifying that this PR works or other notes for the reviewer.

JacksonBurns · 2025-05-25T21:50:37Z

Cantera 2.6 isn't available for Python 3.12, so this PR will also need to upgrade the Cantera version to 3 as mostly completed in #2751

JacksonBurns · 2025-05-26T00:07:45Z

A note for the path forward on this PR - the get_relevant_cycles and get_smallest_set_of_smallest_rings functions will need to be moved to Molecule. Currently they are in Graph, and Graph has no way to send itself to RDKit in order to be subject to RDKit GetSymmSSSR but Molecule can be converted to RDKit via converted.to_rdkit_mol (also exposed as Molecule.to_rdkit_mol, and then RDKit can be used. This may have implications for inheritance elsewhere in the codebase, but should be OK.

jonwzheng · 2025-07-07T23:39:52Z

I ended up moving all the functions in Graph that call get_relevant_cycles or get_smallest_set_of_smallest_rings to Molecule as well. Still have a couple of unit tests to update, but if this is the direction we want to take, hopefully this should be adequate

jonwzheng · 2025-07-16T21:39:51Z

Addressed some of the failures, but still a few nontrivial challenges:

TODO list:

address e / electron representation if integrating to the new RDKit workflow
fix missing vdW bond order
make sure ring algorithms are implemented correctly
certain bond types like R need some sort of RDKit definition; or just never parse them into an RDKit mol object in the first places
how to represent certain surface bonds?
bond orders without definitions? issue from rdkit side or something else?
fill in remaining placeholder unit tests

The tests are just failing mainly for the aromatic compounds, indicating there's some issue with the implementation of the code. Since everything runs, the code "flow" seems clean, just need to iron out the implementation.

jonwzheng · 2025-07-18T19:15:09Z

We are now very close to passing all of the unit tests. There are just a couple major bugs to resolve, namely:

~~how to handle the R/L cutting fragment labels in the RDKit mol object? need to replace with placeholder atom? also with surfaces?~~
~~some bond orders are weird now, so need to investigate where something is going wrong. Due to sanitization maybe?~~
~~I think we can undo the changes to AtomTypeTest~~
~~some surface glitches~~
~~can't kekulize some molecules: if this occurs, maybe back-up is to use increasingly strict flags? e.g. turn kekulization off, and finally turn all sanitize off?~~

jonwzheng · 2025-07-23T21:00:14Z

One of the last failing tests is test_draw_hydrogen_bond_adsorbate, which actually fails due to an assertion error in the codebase:

https://github.com/ReactionMechanismGenerator/RMG-Py/blame/d0d95b3a6bb9cd5d873c8d5b552ec112d6ba9283/rmgpy/molecule/draw.py#L541

@rwest, since you're the last remaining dev who was involved in some capacity with this part of the code, do you know why we have that assert statement? I would think it's better to just plot a "dot" than to outright crash out.
I'm also not sure why this only popped up in this PR - maybe better to address any underlying problem in this code - but thought I'd ask, since this assert might cause problems in the future. There's no safeguards against passing a molecule with no backbone, unless there's something about how this is called that I'm not getting.

While we're at it, the _generate_coordinates also didn't use RDKit to draw if the charge != 0, but maybe that can be relaxed with new versions?

rwest · 2025-07-24T01:43:02Z

One of the last failing tests is test_draw_hydrogen_bond_adsorbate, which actually fails due to an assertion error in the codebase:

https://github.com/ReactionMechanismGenerator/RMG-Py/blame/d0d95b3a6bb9cd5d873c8d5b552ec112d6ba9283/rmgpy/molecule/draw.py#L541

@rwest, since you're the last remaining dev who was involved in some capacity with this part of the code, do you know why we have that assert statement?

I don't. I can see that if there is no straight chain backbone that a function to _find_straight_chain_backbone shouldn't work. It should either fail, and whatever called it deals with the exception (maybe a ValueError is better than AssertError, but either way it's an error to be handled). Or it could return None, but then the calling code would also have to be modified to handle that without crashing. Or the calling code shouldn't try calling it - but then it'd have to know there is no straight chain backbone (and perhaps the easiest way to find that out is to ask for it and fail.) So I think either way the calling code needs updating.

There's no safeguards against passing a molecule with no backbone, unless there's something about how this is called that I'm not getting.

I haven't investigated the full stack trace. But some safeguard, or exception handling, sounds appropriate.

While we're at it, the _generate_coordinates also didn't use RDKit to draw if the charge != 0, but maybe that can be relaxed with new versions?

We can try it. I expect RDKit at the time wasn't being helpful. Hopefully the commit messages offer clues? (this is why we should write helpful commit messages that explain why things are being done).

rwest · 2025-07-24T12:28:26Z

Is this related?

#2744

jonwzheng · 2025-07-24T18:15:26Z

regarding the draw unit test

@rwest yes, thanks, #2744 looks relevant. The test that's failing is for this charged species:

maybe failing due to the missing X-O connectivity, and so no backbone?
I see the discussion in the thread. I'll move any further discussion about this to that issue.
Maybe a separate PR can fix that issue, merged in, and then this can be rebased onto main.

Edit: I tried modifying the connectivity s.t. there is an apparent backbone. But, I still run into this problem even after changing the adjacency list to

1  O u0 p3 c-1 {2,S} {10,H}
2  N u0 p1 c0 {1,S} {3,S} {4,S}
3  O u0 p2 c0 {2,S} {11,S}
4  O u0 p2 c0 {2,S} {7,S}
5  N u0 p1 c0 {6,S} {8,S} {9,S} {7,H}
6  O u0 p2 c0 {5,S} {10,S}
7  H u0 p0 c0 {4,S} {5,H}
8  H u0 p0 c0 {5,S}
9  H u0 p0 c0 {5,S}
10 H u0 p0 c0 {6,S} {1,H}
11 X u0 p0 c0 {3,S}

So it is more fundamentally a problem with the species , than a unit-test specific one.

Note about sanitization

For posterity: some molecules that fail the first kekulization step include:

Can't kekulize mol.  Unkekulized atoms
[H]C1=C([H])C([H])([H])C([H])([H])[c]([H])c1[c]([H])[H]
[H]C([H])=C([H])C1([H])[c]([H])c([H])[c]([H])C1([H])[H]
[H]C([H])=C1c([c]([H])[H])[c]([H])C([H])([H])C1([H])[H]
[H]C1=C([H])C([H])([H])C2([H])[c]([H])c([H])[c]([H])C2([H])C([H])=C1[H]
[H]C1~[C]([H])C2([H])[c](c([H])[c]~1[H])~C([H])~[C]([H])C2([H])C([H])([H])[H]
[H][c]1c([H])[c]([H])C([H])(C([H])([H])[H])C1([H])[H]
[H][c]([H])c1[c]([H])C([H])([H])c2c([H])c([H])c([H])c([H])c-12
[H][C]1~C([H])~[C]([H])C2([H])[C](C([H])([H])[c]([H])c2[c]([H])[H])~C~1[H]
[H][c]([H])c1[c]([H])C([H])([H])c2c([H])c([H])c([H])c([H])c-12
[H][c]1c([H])[c]([H])C([H])(C([H])([H])[H])C1([H])[H]
[H]C1~[C]([H])C2([H])[c](c([H])[c]~1[H])~C([H])~[C]([H])C2([H])[H]
[H][c]([H])c1[c]([H])C1([H])[H]
non-ring atom marked aromatic:
[H][C]~c[c]([H])C([H])([H])C1=C([H])C([H])=C([H])C([H])=C1[H]
[H]c1=c([H])-c([H])=C(C([H])([H])[c]([H])c([H])[c]([H])[H])c([H])=c-1[H]

As implemented in this PR, the RDKit sanitization step will fallback to not use kekulization.

jonwzheng · 2025-07-24T18:42:21Z

Regarding the draw issue, I think I've figured out what's happening...

The test case is drawn as a cyclic molecule, even though it's not truly "cyclic". So it should take the find_cyclic_backbone branch rather than the find_straight_chain_backbone branch.

However, it thinks that len(self.cycles) = 0 because of the changes to ring perception, specifically molecule.get_smallest_set_of_smallest_rings: the H bonds are treated as RDKit hydrogen bonds which don't "count" toward the ring. So in RDKit world it's treated like a non-cyclic species, whereas RMG treats it as cyclic.

So, how to proceed? a few options...

Prefer RDKit drawing: if ring perception is through RDKit, then drawing with RDKit is going to generally be safer. Probably the easiest option, as all we need to do is get rid of the check for if there's charged species. (I think it should work alright personally!)
Change ring perception: instead of treating the H bonds as hydrogen bonds, treat them as an explicit bond?
Create a separate function from SSSR that is used purely to see if there's "ring-like" structure regardless of the nature of the bonds, purely for drawing/graph structure type analysis.

rwest · 2025-07-24T19:52:37Z

I think using RDKit for drawing as much as possible is fine, if it does a good job. Probably we ran into issues in the past that may not persist with more recent versions of RDKit, so it's fine to revisit past decisions. We should continue to put reasons for things in commit messages, to make life easier for future developers (our future selves included).

jonwzheng · 2025-07-24T20:02:26Z

OK. Maybe we can incorporate the fix into #2838 and then rebase onto main once it's merged in.

In ReactionMechanismGenerator#2744 and ReactionMechanismGenerator#2796 it was found that charge-separated bidentate species can have issues due to ring perception conflicts. The previous implementation also by default did not use the rdkit backend for charged species, but this was decided many years ago (~10 years!) In the meantime, RDKit conformer generation has improved and likely this we can just use RDKit by default, which would avoid the pesky edge-case issues for ions/zwitterions. In case the old behavior is desired, use_rdkit can be set to False.

In #2744 and #2796 it was found that charge-separated bidentate species can have issues due to ring perception conflicts. The previous implementation also by default did not use the rdkit backend for charged species, but this was decided many years ago (~10 years!) In the meantime, RDKit conformer generation has improved and likely this we can just use RDKit by default, which would avoid the pesky edge-case issues for ions/zwitterions. In case the old behavior is desired, use_rdkit can be set to False.

JacksonBurns · 2025-08-04T01:52:40Z

With the merger of e92bec0 this PR's CI should now pass - let's see!

JacksonBurns · 2025-08-04T13:30:14Z

@jonwzheng the overnight test failures look spurious - re-running them now.

If you had remove_h = True and sanitize="partial", it would try passing that along to RDKit, which would not expect it. Now we remove H's (without sanitizing) before possibly doing the sanitizing.

So you don't need to create an instance of Fragment() just in order to call the method (seeing as that seems to be how it is mostly used.) Probably it should be a module level function, but this is a minimal change.

Now that we use RDKit to detect rings, we create a lot of rdkit molecules in which we don't care about the bond orders, sometimes we have non-integer or weird bond orders, and we also don't want to waste time sanitizing them since all we care about is connectivity of the molecular graph. So if you know you're calling to_rdkit_mol in this setting, call it with ignore_bond_orders=True and sanitize=False.

This should save complication and time. Also, don't sanitize. All we're going to ask it for is ring information. For that all bonds (that we give it) count. If we want it to ignore H-bonds, or vdW bonds, for some scenarios, then we should figure that out. But often we use this ring detection code for things like drawing, which does want to use those bonds.

This reverts commit 55f8fc4. We should, during ring detection, now be not attempting to convert bond orders into strings, so this warning should not be triggered. We can restore the debug level if this becomes too noisy.

We do a call to detect_cutting_label (which is a complicated function to find them buried inside SMILES strings) passing it the label of every atom. (Though maybe these are mostly ''?). This commit checks whether a simpler approach would work: just checking if the atom.label is either 'R' or 'L'. If that's the case, we can then make the change.

As far as I can tell, cutting labels can only be "L" or "R". I had a previous commit that checked this is equivalent, and by running the unit tests and regression tests could not find a counterexample. This check is much simpler, and about 200 times faster, than the previous call to detect_cutting_label. Also simplified the dictionary structure.

…rnings. Previously, if you call to_rdkit_mol on a Fragment, it would always return the mapping, but if you called it on a Molecule, it would by default not. Now if you call it with return_mapping=False (or not specified) then it doesn't return the mapping, and if you call it with return_mapping=True, then it does. Internally, it calls converter.to_rdkit_mol with return_mapping=True, but doesn't pass it back unless asked. Now the behavior of thing.to_rdkit_mol() is the same whether thing is a Molecule or a Fragment (which is a subclass of Molecule). Should reduce number of warnings, eg. when drawing fragments. Also cleaned up the code a bit, and added some docstrings,

…_rdkit_mol The latest change of to_rdkit_mol means that both Molecule and Fragment return the same signature, and the calling code doesn't need to cope with unexpected tuples.

The comment suggests it was going through Geometry in order to save the mapping, but this can be done with a direct call to to_rdkit_mol using return_mapping=True.

- cache the result of the call to rdkit to avoid remaking the molecule and running the decomposition - use GetSSSR only and deprecate other methods (and do so following conventions correctly) - add a test for the OPTIONAL symmetric SSSR passthrough

JacksonBurns · 2026-02-13T16:39:12Z

RMS install fails on 3.10 and 3.11 because of something to do with PythonCall and RMS

JacksonBurns · 2026-02-13T16:49:51Z

RMS install fails on 3.10 and 3.11 because of something to do with PythonCall and RMS

Based on this:

ERROR: Unsatisfiable requirements detected for package PythonCall [6099a3de]: PythonCall [6099a3de] log: ├─possible versions are: 0.1.0-0.9.31 or uninstalled ├─restricted to versions 0.9.31 by project [b7569dc2], leaving only versions: 0.9.31 │ └─project [b7569dc2] log: │ ├─possible versions are: 0.0.0 or uninstalled │ └─project [b7569dc2] is fixed to version 0.0.0 └─restricted to versions 0.9-0.9.28 by ReactionMechanismSimulator [c2d78dd2] — no versions left └─ReactionMechanismSimulator [c2d78dd2] log: ├─possible versions are: 1.0.1 or uninstalled └─ReactionMechanismSimulator [c2d78dd2] is fixed to version 1.0.1

The solution may be to increment Julia version to 1.10 or newer, which would allow us to use the latest versions of PythonCall (that support more versions of Python), and is the minimum supported version according to PythonCall docs.

…nCall version setting

JacksonBurns · 2026-02-13T17:24:37Z

On second thought, Julia actually already is 1.10 - the problem appears to be that when we install pyjuliacall it locks the version of PythonCall. If we install the former after the latter, it may resolve it.

JacksonBurns · 2026-02-13T17:30:10Z

@copilot - the changes in this pull request are mostly complete, except that these two tests are failing:

FAILED test/rmgpy/data/thermoTest.py::TestThermoDatabase::test_thermo_estimation_not_affect_database - AssertionError: assert 'PolycyclicRing' in []
FAILED test/rmgpy/data/thermoTest.py::TestCyclicThermo::test_thermo_for_monocyclic_and_polycyclic_same_molecule - assert 0 == 1
 +  where 0 = len([])

There's something wrong with the detection of polycycles, though the detection of the smallest set of smallest rings appears correct. Can you figure out what is causing this issue, and fix things so that these last two unit tests pass?

rmgpy/molecule/molecule.py

JacksonBurns marked this pull request as draft May 25, 2025 21:43

JacksonBurns force-pushed the feat/rdkit_sssr branch from 99bf00c to 8d6fe64 Compare May 25, 2025 21:49

JacksonBurns mentioned this pull request May 29, 2025

Add T-Dependent and ML Options for Solvation Corrections #2798

Closed

Base automatically changed from feat/py39_rebase to main May 29, 2025 19:55

jonwzheng self-assigned this Jun 30, 2025

jonwzheng mentioned this pull request Jun 30, 2025

Update Jupyter notebooks and associated codebase [RMG-Py] #2815

Open

13 tasks

JacksonBurns mentioned this pull request Jul 9, 2025

Make RMG-Py conda Work Again 🎉 #2820

Merged

jonwzheng mentioned this pull request Jul 24, 2025

Updates to draw.py are causing a crash for charge separated bidentates #2744

Closed

jonwzheng mentioned this pull request Jul 25, 2025

Fix drawing errors for charged species by moving to rdkit geometry generation #2838

Merged

jonwzheng force-pushed the feat/rdkit_sssr branch from 7598c98 to 019362c Compare July 25, 2025 17:09

JacksonBurns force-pushed the feat/rdkit_sssr branch from 6e6490a to e92bec0 Compare August 4, 2025 01:51

JacksonBurns mentioned this pull request Aug 4, 2025

Status Tracking: chemprop and ML Thermochemistry Prediction #2559

Open

1 task

rwest and others added 22 commits February 13, 2026 11:08

Fix sanitization issue in to_rdkit_mol

5c0cbcb

If you had remove_h = True and sanitize="partial", it would try passing that along to RDKit, which would not expect it. Now we remove H's (without sanitizing) before possibly doing the sanitizing.

Make detect_cutting_label a static method.

3ee940c

So you don't need to create an instance of Fragment() just in order to call the method (seeing as that seems to be how it is mostly used.) Probably it should be a module level function, but this is a minimal change.

Revert "Change debug level of RDKit-related warnings"

99007e6

This reverts commit 55f8fc4. We should, during ring detection, now be not attempting to convert bond orders into strings, so this warning should not be triggered. We can restore the debug level if this becomes too noisy.

Ring finding code doesn't need to cope with unwanted mappings from to…

d01b57f

…_rdkit_mol The latest change of to_rdkit_mol means that both Molecule and Fragment return the same signature, and the calling code doesn't need to cope with unexpected tuples.

Tweak to MoleculeDrawer: don't bother making a Geometry object.

3461488

The comment suggests it was going through Geometry in order to save the mapping, but this can be done with a direct call to to_rdkit_mol using return_mapping=True.

use sssr instead of symmsssr

e66d0e0

disable build isolation for compat with py 3.10 and 3.11

327bc68

explicitly disable build isolation

bfe48a3

upper bound python version

684fcbb

the mac0s-13 runner no longer exists

af20498

refactor to use getsssr, see extended

01595f0

- cache the result of the call to rdkit to avoid remaking the molecule and running the decomposition - use GetSSSR only and deprecate other methods (and do so following conventions correctly) - add a test for the OPTIONAL symmetric SSSR passthrough

update header

5636a4b

reminder to self: engage brain before adding attributes

19025f9

change attr setting

f75a635

update header for cache so attribute is allocated

31a827b

fix pytest syntax

ca71778

update expected error type

1aa6de5

JacksonBurns force-pushed the feat/rdkit_sssr branch from 68fbf96 to 1aa6de5 Compare February 13, 2026 16:10

install pyjuliacall after installing rms, which allows flexible Pytho…

af8f78d

…nCall version setting

rwest reviewed Feb 13, 2026

View reviewed changes

rmgpy/molecule/molecule.py Show resolved Hide resolved

setup python side of RMS after complete rms install done

c5ac7aa

Conversation

JacksonBurns commented May 25, 2025

Motivation or Problem

Description of Changes

Testing

Reviewer Tips

Uh oh!

JacksonBurns commented May 25, 2025

Uh oh!

JacksonBurns commented May 26, 2025

Uh oh!

jonwzheng commented Jul 7, 2025

Uh oh!

jonwzheng commented Jul 16, 2025 • edited by JacksonBurns Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jonwzheng commented Jul 18, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jonwzheng commented Jul 23, 2025

Uh oh!

rwest commented Jul 24, 2025

Uh oh!

rwest commented Jul 24, 2025

Uh oh!

jonwzheng commented Jul 24, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

regarding the draw unit test

Note about sanitization

Uh oh!

jonwzheng commented Jul 24, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

rwest commented Jul 24, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jonwzheng commented Jul 24, 2025

Uh oh!

JacksonBurns commented Aug 4, 2025

Uh oh!

JacksonBurns commented Aug 4, 2025

Uh oh!

JacksonBurns commented Feb 13, 2026

Uh oh!

JacksonBurns commented Feb 13, 2026

Uh oh!

JacksonBurns commented Feb 13, 2026

Uh oh!

JacksonBurns commented Feb 13, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

jonwzheng commented Jul 16, 2025 •

edited by JacksonBurns

Loading

jonwzheng commented Jul 18, 2025 •

edited

Loading

jonwzheng commented Jul 24, 2025 •

edited

Loading

jonwzheng commented Jul 24, 2025 •

edited

Loading

rwest commented Jul 24, 2025 •

edited

Loading