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

[kirkbadger18]

When I try to run RMG in cases where charge separated bidentates show up, the recent updates to draw.py cause the following error:

I have isolated the issue by running the following input file in rmg-py main and rmg-database main with just the addition of one reaction in a reaction library. This one reaction contains a charge separated bidentate: XOXNO.

The dictionary looks like:

And the reaction looks like:

Here are the files needed to regenerate this error:
generate_error.zip

I think it has to do with this check in draw.py:

I tried deleting this check and it fixed the error for me, but I'm sure these lines must be here for a reason, and maybe deleting them would cause issues for someone else.

[kirkbadger18]

@rwest, this is the draw.py issue we discussed.

[rwest]

On a new branch, please can you add a (failing) unit test in test/rmgpy/molecule/drawTest.py
(eg. like test_draw_bidentate_adsorbate )
and open a pull request from that branch?
That will help us figure out what is crashing _find_straight_chain_backbone.

I wasn't clear from your description what your modification is doing, but from the stack trace it looks like the problem was with using our own coordinate generation, so your change is probably making it use RDKit instead.
It looks like the test was previously saying if any of the atoms have a charge then don't use RDKit.
I think introduced in ad80018 with the note "Draw molecule using own algorithm if one of the molecule's atoms is charged. This has more likelihood of not failing on rdkit due to incorrect valence."
Perhaps things have changed in RDKit in the intervening 10 years (!!) and it doesn't complain or fail any more.
Your fix may be fine, but to be safe, we should ensure there are lots of test cases for drawing zwitterions and perhaps even ions. So maybe also add some of these to the drawTest.py?

[github-actions]

This issue is being automatically marked as stale because it has not received any interaction in the last 90 days. Please leave a comment if this is still a relevant issue, otherwise it will automatically be closed in 30 days.

[rwest]

@kirkbadger18 should this be reopened and not stale?

[kirkbadger18]

This is still an issue for me, I just forgot about it since my fix worked for me.

[kirkbadger18]

I'll take a shot at writing this test, I think I can just copy the format here

[kirkbadger18]

@rwest
I just opened #2838 , let's see if it passes all its tests but the one I added first, the linter always seems to get me.

[jonwzheng]

I personally vote for using rdkit to do drawing by default for molecules, and keeping use_rdkit as an argument.

in #2796, we're updating the molecule code to use RDKit for ring perception rather than ringdecomposerlib on RMG molecule objects. I encountered a (similar) crash to this case because RDKit inferred that the species was non-cyclic, but the RMG-based draw algorithm thinks it's cyclic in draw.

I outlined a few possible solutions in the thread there:

1. 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! I have had pretty good experiences with using RDKit for charged species. As Richard said we probably would want to incorporate some unit tests just to be safe, and maybe also manual examination.)
2. Change ring perception: instead of treating the H bonds as hydrogen bonds, treat them as an explicit bond? (This probably will lead to unintended behavior since we don't really want to perceive H-bonded complexes as "rings")
3. 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.

[jonwzheng]

BTW, the problem specifically for your test case @kirkbadger18 is that during drawing, the molecule becomes cyclic (because the surface sites are connected, _connect_surface_sites before generate_coordinates()). As far as I can tell, this is done just to make the drawn result look better.

But, because ring perception is performed before the surface sites are connected, the molecule is inferred to be non-cyclic in generate_coordinates. And therefore it tries to find a backbone, but it fails because there's only one terminal site (the O sticking out; everything else is part of the "artificial" cycle).

To fix this, we would need to refactor the drawing logic. I think during connect_surface_sites you could give the surface sites a special order bond (e.g. like 0.1 or 0.5 or something). Then, in _generate_straight_chain_coordinates, check for terminal atoms only checking among bonds with bond order 1.0.

Or, alternatively, we can forget it and just switch everything to the RDKit backend, and make a note that we don't maintain the non-RDKit method.

[kirkbadger18]

@jonwzheng Thanks for looking into this! I do not have a preference on how it should be solved; I just wanted to mention that this was happening. @rwest what are your thoughts?

[JacksonBurns]

+1 for deprecating our internal method and switching to rdkit entirely