Isomorphism comparisons are slow on similar large molecules

[galioguo]

Description and Output

I have two quite large (291 atoms each) stereoisomers that differ at precisely one chiral site. I'll call them offmol1 and offmol2. I want to feed them into an Interchange using from_smirnoff. I precomputed their partial charges, and fed this in through the 'charge_from_molecules' field.

import logging
import openff.toolkit as tk

logging.basicConfig(level=logging.DEBUG)
topology = tk.Topology.from_molecules([offmol1, offmol2])
sage = tk.ForceField('openff_unconstrained-2.2.1.offxml')
interchange = Interchange.from_smirnoff(
    force_field = sage, 
    topology = topology, 
    charge_from_molecules = [offmol1, offmol2], 
    allow_nonintegral_charges = True
)

This code hangs indefinitely. Upon checking debug logs, this originates from the isomorphism checks under _assign_charges_from_molecules inside interchange.smirnoff._nonbonded.py

def _assign_charges_from_molecules(
    cls,
    topology: Topology,
    unique_molecule: Molecule,
    charge_from_molecules=Optional[list[Molecule]],
) -> tuple[bool, dict, dict]:
    if charge_from_molecules is None:
        return False, dict(), dict()

    for molecule_with_charges in charge_from_molecules:
        if molecule_with_charges.is_isomorphic_with(unique_molecule):
            break
    else:
        return False, dict(), dict()

    _, atom_map = Molecule.are_isomorphic(
        molecule_with_charges,
        unique_molecule,
        return_atom_map=True,
    )

Although partial charges are in general dependent on stereochemistry, it would be helpful to add an option to ignore stereochemistry when specifying charges from molecules.

[mattwthompson]

Only have a minute to look at this now, will be able to spend more time on it tomorrow - does something like offmol1.is_isomorphic_with(offmol2) also hang indefinitely?

Also if you're able to share files/SMILES/etc to reproduce this that would be helpful, I would like to be able to run it on my machine

[galioguo]

Yes, offmol1.is_isomorphic_with(offmol2) seems to hang indefinitely as well. Looks like it's stuck in a a cycle of recursive calls inside NetworkX.

I attached two TXT files (Github won't let me upload MOL files, and I don't think the Topology.from_molecules can eat TXT format, so it might be necessary to change the extension to MOL) and some code that reproduces the issue. I appreciate your prompt response.

import openff.toolkit as tk
from openff.interchange import Interchange
import logging

logging.basicConfig(level=logging.DEBUG)
offmol1 = tk.Molecule.from_file('./offmol1.mol')
offmol2 = tk.Molecule.from_file('./offmol2.mol')
offmol1.assign_partial_charges('gasteiger')
offmol2.assign_partial_charges('gasteiger')

topology = tk.Topology.from_molecules([offmol1, offmol2])
sage = tk.ForceField('openff_unconstrained-2.2.1.offxml')
interchange = Interchange.from_smirnoff(
    force_field = sage, 
    topology = topology, 
    charge_from_molecules = [offmol1, offmol2], 
    allow_nonintegral_charges = True
)

offmol1.txt
offmol2.txt

[mattwthompson]

This turns out to be a very messy problem! You correctly diagnosed the fulcrum of the slowdown. Unfortunately, I can't come up with a safe and general solution right now.

• The toolkit is slow to run isomorphism checks of large molecules
• Your molecules differ in only one bit and the difference is somewhere in the middle of the molecule, which is particularly hard for the toolkit to deal with
• Isomorphism checks on pre-set molecules can't safely be skipped in general; even though in your case the stereochemistry doesn't impact partial charges, it often does
• Exposing options for isomorphism checks in pre-set molecules is a scale of complexity that I'm not eager to add to the public API. You may guess that a couple of arguments passed through at the right place will get this working, and you'd be right, but that's one use case in a large array of options that would need to be supported.

Are there any other molecules in your system? I can put together a hacky workaround that would apply Gastiger charges to all molecules, except for water and monoatomic ions.

[galioguo]

Thanks for your reply! This does seem like a messy problem to fix in full generality.

Since a particular quirk of my system is that chiral centers mostly have a hydrogen present, I made the generally inaccurate approximation that partial charges don't depend on stereochemistry. My engineering fix was to ignore stereochemistry during isomorphism checks, making the following changes in my installation of OpenFF:

~/galio/.conda/envs/polymers9/lib/python3.10/site-packages/openff/interchange/smirnoff/_nonbonded.py
    LINES 791, 801, added additional argument, 'atom_stereochemistry_matching = False'

~/galio/.conda/envs/polymers9/lib/python3.10/site-packages/openff/toolkit/topology/topology.py
    LINES 1198, added additional argument, 'atom_stereochemistry_matching = False'

Although this is somewhat ad-hoc, it allows me to construct the interchange.

My system contains other molecules that are much larger (~1000 atoms) and have similar chemistry (I am simulating a crosslinked polymer). My hope is that 'charge_from_molecules' will be robust enough to handle even larger molecules. I would definitely be interested in seeing what kind of workaround you're thinking of. Would your approach allow for applying charges from other charge methods?

[mattwthompson]

I'm going to move this into the toolkit repo since the problem is with slow isomorphism checks and there isn't much Interchange can do to work around it in general. #353 is a good stub issue and the polymer-performance label also has some similar issues.

I'll introduce my hacky idea in a separate issue within the week (please ping me if I don't get around this) just to keep this one focused