Cenverting from RDKit molcule to OpenMM leads to incompatibility

[MiloszGrabski]

Describe the bug
When adding a molecule converted from RDKit using OpenFF Toolkit to an OpenMM Modeller, the following error is raised:

AttributeError: 'numpy.ndarray' object has no attribute 'x'

This occurs due to a mismatch in coordinate format returned by t.get_positions()—a NumPy array— instead Vec3

Minimal example

This helper function causes the issue:

def add_molecule_to_modeller(modeller: Modeller, molecule: rdkit.Chem.Mol, resname: str = "FRA") -> None:
    resid = modeller.topology.getNumResidues() + 1
    mol = Molecule.from_rdkit(molecule)
    t = mol.to_topology()
    for atom in t.atoms:
        atom.metadata["residue_name"] = resname
        atom.metadata["residue_number"] = str(resid)
    modeller.add(t.to_openmm(), t.get_positions().to_openmm())  # <-- Causes error

Downstream, this triggers an exception when trying to compute bounding box coordinates after using modeller.addSolvent()

    for pos in water_positions:
        min_coords = Vec3(min(min_coords.x, pos.x), min(min_coords.y, pos.y), min(min_coords.z, pos.z))
        max_coords = Vec3(max(max_coords.x, pos.x), max(max_coords.y, pos.y), max(max_coords.z, pos.z))

Expected behavior
The to_openmm() call should produce a list of Vec3 positions compatible with Modeller.addSolvent().

Actual behavior
The t.get_positions().to_openmm() call returns a NumPy array, not a list of Vec3, leading to the AttributeError.

Workaround
Manually convert coordinates to Vec3:

from openmm.unit import nanometer, Vec3

def add_molecule_to_modeller(modeller: Modeller, molecule: rdkit.Chem.Mol, resname: str = "FRA") -> None:
    """Converts RDKit molecule to topology and adds it to the modeller."""
    resid = modeller.topology.getNumResidues() + 1
    mol = Molecule.from_rdkit(molecule)
    t = mol.to_topology()
    for atom in t.atoms:
        atom.metadata["residue_name"] = resname
        atom.metadata["residue_number"] = str(resid)
    positions = []
    for x, y, z in t.get_positions():
        x = x.to('nanometer').magnitude
        y = y.to('nanometer').magnitude
        z = z.to('nanometer').magnitude
        positions.append(Vec3(x, y, z) * nanometer)
    modeller.add(t.to_openmm(), positions)

Computing environment:

openff-toolkit 0.16.9 pyhd8ed1ab_0 conda-forge
openmm 8.2.0 py312haeaa2e4_2 conda-forge

[mattwthompson]

For better or worse, openff.toolkit.Quantity.to_openmm() returns an openmm.unit.Quantity which itself attempts to wrap the raw numbers in a similar way to how they're wrapped with OpenFF's unit solution. Always converting to Vec3 is a non-starter because quantities can also represent things like force constants, constraint distances, etc.

Your chaining together of these two methods make this result seem strange, but keep in mind that .get_positions doesn't know whether the positions will be used with OpenMM or some other application later. I think this is just a funky quirk of OpenMM - sometimes openmm.unit.Quantity-wrapped NumPy arrays are okay, sometimes they're not.

That being said, I would be frustrated if I was in your shoes. Dealing with Vec3 when I basically have the right information is a hassle. I wonder if this justifies a sort of Quantity.to_vec3 method? Or possibly a Topology.get_vec3_positions? I would like to bounce this idea off of @Yoshanuikabundi

[Yoshanuikabundi]

Hmm this is definitely bad ergonomics. If I were in charge of everything, I think the way I'd fix this is to fix Modeller so that it can accept an array - I usually expect to be able to use arrays and Vec3s interchangeably with OpenMM tools. You might like to raise this as an issue in OpenMM and see if they'd be interested in a fix? It might just be an oversight.

We could detect if the last dimension of the array has length 3, and then Vec3-ify only in that case... but this would be a really strange behavior and I expect it would break more than it would fix. As @mattwthompson said, the purpose of the Quantity.to_openmm method is to convert an OpenFF quantity to OpenMM, not to generally convert any object wrapped in a quantity to an object that OpenMM is guaranteed to be able to read.

There are methods on the upstream Pint Quantity class that do provide ergonomics for converting between lists of quantities and quantities of lists, quantity sequences to arrays, and so on, so a Quantity.to_vec3 method wouldn't be the craziest thing we could possibly do. But if this is just an OpenMM oversight, it'd be easier to get it fixed there - that way users don't need to learn anything new and the to_openmm method just gets its magic back.

[MiloszGrabski]

Hmm this is definitely bad ergonomics. If I were in charge of everything, I think the way I'd fix this is to fix Modeller so that it can accept an array - I usually expect to be able to use arrays and Vec3s interchangeably with OpenMM tools. You might like to raise this as an issue in OpenMM and see if they'd be interested in a fix? It might just be an oversight.

Modeller is fine, the problem is how those objects get added when using addSolvent(). I only noticed the problem when tried using np.max. Normally, I would not even notice that data types are different as everything works fine. Then you have type mismatch as water position values are different type than ones added using OpenFF.

[j-wags]

Thanks for the report and the good discussion here. Agree that this is obnoxious but I also don't see a single technical fix that would clear this all up without adding negative side effects. One step in the right direction could be adding a code example to the to_openmm docstring showing how to convert the output to vec3s?

[mattwthompson]

IMHO this is squarely a problem of OpenMM's API, but an example bridging this gap wouldn't hurt

[MiloszGrabski]

Either example or comment about potential variable datatypes. I create custom solvent boxes myself, so I am used to OpenMM being highly sensitive to units and data types, but for someone who just started it would be infuriating. Hence, anything that would make people life easier would be of great benefit.