Catch/raise descriptive error on conformer generation failure

[mattwthompson]

Kinda addresses #1049, at least for the RDKit case.

• Tag issue being addressed
• Add tests
• Update docstrings/documentation, if applicable
• Lint codebase
• Update changelog

[codecov]

Codecov Report

Merging #1050 (11d3bdf) into master (7a3fad5) will decrease coverage by 0.33%.
The diff coverage is 100.00%.

[lgtm-com]

This pull request introduces 1 alert when merging 1c79ec1 into f41e671 - view on LGTM.com

new alerts:

• 1 for Unused import

[lilyminium]

Does RDKit output any of its own error messages?

[mattwthompson]

Not as far as I'm aware, but I don't use RDKit directly too often so I could be missing something.

>>> from rdkit import Chem
>>> from rdkit.Chem import AllChem
>>> rdmol = Chem.MolFromSmiles(100 * "CC")
>>> rdmol = Chem.AddHs(rdmol)
>>> AllChem.EmbedMultipleConfs(rdmol, numConfs=1, pruneRmsThresh=1)
<rdkit.rdBase._vecti object at 0x13763e200>
>>> rdmol.GetNumConformers()
0

[lilyminium]

That's disappointing :(.

Maybe you could query the output directly? It's supposed to contain the conformer IDs of the conformers generated, so if there aren't any, that should correspond to conformer generation failing. Conformers are already cleared by default (clearConfs=True) so this should lead to the same behaviour.

status = AllChem.EmbedMultipleConfs(
            rdmol,
            numConfs=n_conformers,
            pruneRmsThresh=rms_cutoff / unit.angstrom,
            randomSeed=1,
            # params=AllChem.ETKDG()
        )
if not status:
    raise ConformerGenerationError("RDKit conformer generation failed.")

Incidentally, I do get a conformer if I start off with random coordinates. I wonder if the distance matrix is having difficulty with the number of atoms.

From the docs:

useRandomCoords: Start the embedding from random coordinates instead of using eigenvalues of the distance matrix.

>>> status = AllChem.EmbedMultipleConfs(rdmol, numConfs=1, useRandomCoords=True)
>>> list(status)
[0]

Looks alright I guess.

[mattwthompson]

Looks like a reasonable C200 to me.

I don't think this is a good way of testing for this failures, since other settings or future RDKit updates could make generating stupidly large molecules safer/faster, and I'd like to avoid asking it to do an even larger molecule like a million-carbon alkane just to catch a failure. I figure there's got to be a small molecule out there that will reliably fail to be turned into 3D, but not for other reasons (i.e. RDKit wouldn't throw an error for other reasons). But that's as far as I got with that line of thinking. Any suggestions?

[mattwthompson]

One option would be to blindly grab something from their issue tracker, but that's more or less the opposite of maintainable code. The C200 I have right now fails in about 50 seconds on my machine, but two of the smaller molecules from rdkit/rdkit#1433 fail in about 15 seconds each.

[lilyminium]

No idea, sorry 😬 from that issue I'd guess it's flexible or large molecules that give it trouble, although I can't work out what's especially hard about "C1CC1c2[nH]nc3C4CCC4c23". I did make a variation on a theme ("C12C3C(C1C2)C4=N[N](C(=C34)C5C6C5CC6)C7CC7") for something that might not get directly addressed in the future.

>>> rdmol = Chem.MolFromSmiles("C12C3C(C1C2)C4=N[N](C(=C34)C5C6C5CC6)C7CC7")
>>> rdmol = Chem.AddHs(rdmol)
>>> status = AllChem.EmbedMultipleConfs(rdmol, numConfs=1, useRandomCoords=True)
>>> len(status)
0

[mattwthompson]

I asked in Slack and the PIs chimed in with a list of suggestions, the smallest of which was uranium hexafluoride, which I'm rolling with. Coincidentally it has the funniest SMILIES string I've seen to date ...

[mattwthompson]

While working on the unit tests, I came across some molecules that OEChem doesn't like but RDKit likes (enough to generate something):

>>> from openff.toolkit.topology import Molecule
>>> smi = "B1(OB(OB(O1)OC)OC)OC"
>>> molecule = Molecule.from_smiles(smi)
>>> from openff.toolkit.utils.toolkits import OpenEyeToolkitWrapper, RDKitToolkitWrapper, AmberToolsToolkitWrapper, GLOBAL_TOOLKIT_REGISTRY, ToolkitRegistry
>>> molecule.generate_conformers(n_conformers=1, toolkit_registry=RDKitToolkitWrapper())
>>> molecule.generate_conformers(n_conformers=1, toolkit_registry=OpenEyeToolkitWrapper())
Warning: no MMFF bend parameters for atoms 1 2 3
Warning: Unable to find bend parameters for
Warning: no MMFF bend parameters for atoms 1 2 3
Warning: : Force field setup failed due to missing parameters
Warning: no MMFF bend parameters for atoms 1 2 3
Warning: Unable to find bend parameters for
Warning: no MMFF bend parameters for atoms 1 2 3
Warning: : Force field setup failed due to missing parameters
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/mwt/software/openforcefield/openff/toolkit/topology/molecule.py", line 3230, in generate_conformers
    return toolkit.generate_conformers(
  File "/Users/mwt/software/openforcefield/openff/toolkit/utils/openeye_wrapper.py", line 1693, in generate_conformers
    raise ConformerGenerationError(
openff.toolkit.utils.exceptions.ConformerGenerationError: OpenEye Omega conformer generation failed
>>> molecule.generate_conformers(n_conformers=1)
Warning: no MMFF bend parameters for atoms 1 2 3
Warning: Unable to find bend parameters for
Warning: no MMFF bend parameters for atoms 1 2 3
Warning: : Force field setup failed due to missing parameters
Warning: no MMFF bend parameters for atoms 1 2 3
Warning: Unable to find bend parameters for
Warning: no MMFF bend parameters for atoms 1 2 3
Warning: : Force field setup failed due to missing parameters
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/mwt/software/openforcefield/openff/toolkit/topology/molecule.py", line 3221, in generate_conformers
    return toolkit_registry.call(
  File "/Users/mwt/software/openforcefield/openff/toolkit/utils/toolkit_registry.py", line 366, in call
    raise e
  File "/Users/mwt/software/openforcefield/openff/toolkit/utils/toolkit_registry.py", line 362, in call
    return method(*args, **kwargs)
  File "/Users/mwt/software/openforcefield/openff/toolkit/utils/openeye_wrapper.py", line 1693, in generate_conformers
    raise ConformerGenerationError(
openff.toolkit.utils.exceptions.ConformerGenerationError: OpenEye Omega conformer generation failed
GLOBAL_TOOLKIT_REGISTRY
ToolkitRegistry containing OpenEye Toolkit, The RDKit, AmberTools, Built-in Toolkit
>>> molecule.generate_conformers(n_conformers=1, toolkit_registry=GLOBAL_TOOLKIT_REGISTRY)
Warning: no MMFF bend parameters for atoms 1 2 3
Warning: Unable to find bend parameters for
Warning: no MMFF bend parameters for atoms 1 2 3
Warning: : Force field setup failed due to missing parameters
Warning: no MMFF bend parameters for atoms 1 2 3
Warning: Unable to find bend parameters for
Warning: no MMFF bend parameters for atoms 1 2 3
Warning: : Force field setup failed due to missing parameters
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/mwt/software/openforcefield/openff/toolkit/topology/molecule.py", line 3221, in generate_conformers
    return toolkit_registry.call(
  File "/Users/mwt/software/openforcefield/openff/toolkit/utils/toolkit_registry.py", line 366, in call
    raise e
  File "/Users/mwt/software/openforcefield/openff/toolkit/utils/toolkit_registry.py", line 362, in call
    return method(*args, **kwargs)
  File "/Users/mwt/software/openforcefield/openff/toolkit/utils/openeye_wrapper.py", line 1693, in generate_conformers
    raise ConformerGenerationError(
openff.toolkit.utils.exceptions.ConformerGenerationError: OpenEye Omega conformer generation failed

As far as I understand it, the first two calls do what we expect; tell the toolkit which wrapped toolkit to use, and report those results. However, I expected the vanilla case to fail on OpenEye but then proceed to try RDKit. This is what .assign_partial_charges does, or so I thought, and I don't see why it's not behaving the same way.

[lilyminium]

You would need to set raise_exception_types=[], (or in some way exclude the ConformerGenerationError), in molecule.generate_conformers for the fall-through, like assign_partial_charges does.

[lilyminium]

It might have been easier to define an allow-list of exceptions, really...

[mattwthompson]

Ah, you're right. I forgot the falsiness of the default value

openff-toolkit/openff/toolkit/utils/toolkit_registry.py

Lines 354 to 355 in f41e671

	if raise_exception_types is None:
	raise_exception_types = [Exception]

I agree there might be a better way to control this behavior, this behavior resulting from passing an empty list is not intuitive to me.

Behaving as expected now, I think:

In [1]: >>> from openff.toolkit.topology import Molecule
   ...: >>> smi = "B1(OB(OB(O1)OC)OC)OC"
   ...: >>> molecule = Molecule.from_smiles(smi)
   ...: >>> from openff.toolkit.utils.toolkits import OpenEyeToolkitWrapper, RDKitToolkitWrapper, AmberToolsToolkitWrapper, GLOBAL_TOOLKIT_REGISTRY, ToolkitRegistry
   ...:

In [2]: >>> molecule.generate_conformers(n_conformers=1)
   ...:
Warning: no MMFF bend parameters for atoms 1 2 3
Warning: Unable to find bend parameters for
Warning: no MMFF bend parameters for atoms 1 2 3
Warning: : Force field setup failed due to missing parameters
Warning: no MMFF bend parameters for atoms 1 2 3
Warning: Unable to find bend parameters for
Warning: no MMFF bend parameters for atoms 1 2 3
Warning: : Force field setup failed due to missing parameters

[j-wags]

Wow, thanks @lilyminium and @mattwthompson for digging into this! Yeah, the exception blocklist was a weird call I made during an attempt to appease the ToolkitAM1BCCHandler as I was trying to teach it what is vs. isn't fatal during system creation. This PR is a massive improvement to a lot of things, and past some non-blocking comments, the only additional thing I'd really like to see is a releasenotes update for the behavior change of molecule.generate_conformers.