openforcefield/openff-toolkit

Allow `Topology.set_positions(None)`

mattwthompson opened this issue · 3 comments

mattwthompson commented

Is your feature request related to a problem? Please describe.

I'm dealing with a use case in which I want to modify a topology's positions but reset them to what they were before. The getter and setter introduced in #1368 makes this clean and first-class except for the (common, I'd argue) case in which there aren't originally positions at all.

Here is an example of how this fails and what it takes to work around it:

import numpy
from openff.toolkit import Molecule, Quantity, unit
from openff.toolkit.utils.exceptions import IncompatibleUnitError


def temporarily_modify(topology, new_positions):
    """Modify a topology's positions, but then reset them."""
    original_positions = topology.get_positions()
    topology.set_positions(new_positions)
    print(f"Modified positions mean: {numpy.mean(topology.get_positions().m)}")
    topology.set_positions(original_positions)
    print(f"Reset positions mean: {numpy.mean(topology.get_positions().m)}")


topology = Molecule.from_smiles("CCO").to_topology()

ones = Quantity(
    numpy.ones(shape=(topology.n_atoms, 3)),
    unit.nanometer,
)

try:
    temporarily_modify(topology, 3 * ones)
except IncompatibleUnitError as error:
    # Fails because the original positions (None) can't be set
    print("Can't reset positions: " + str(error))
    # Prints
    # Can't reset positions: array should be an OpenFF Quantity with dimensions of length

# Can work around this by only calling it on a topology with positions
molecule = Molecule.from_smiles("CCO")
molecule._conformers = [0 * ones]
topology_with_positions = molecule.to_topology()

# Works; prints 4.0 and 0.0
temporarily_modify(topology_with_positions, 4.0 * ones)

# or can work around by special-casing downstream, which is not desirable


def workaround(topology, new_positions):
    """Modify a topology's positions, but then reset them."""
    original_positions = topology.get_positions()
    topology.set_positions(new_positions)
    print(f"Modified positions mean: {numpy.mean(topology.get_positions().m)}")
    if original_positions is None:
        for molecule in topology.molecules:
            molecule._conformers = None
        print(f"Reset positions mean: {topology.get_positions()}")
    else:
        topology.set_positions(original_positions)
        print(f"Reset positions mean: {numpy.mean(topology.get_positions().m)}")


topology = Molecule.from_smiles("CCO").to_topology()

# Works; prints 5.0 and 0.0
workaround(topology, 5.0 * ones)

Running it gives the output:

$ python modify.py
Modified positions mean: 3.0
Can't reset positions: array should be an OpenFF Quantity with dimensions of length
Modified positions mean: 4.0
Reset positions mean: 0.0
Modified positions mean: 5.0
Reset positions mean: None

Describe the solution you'd like

Just let array=None, which is processed internally along the lines of [molecule._conformers = None for molecule in self.molecules].

Describe alternatives you've considered

See above and below.

Additional context

I think it's reasonable for users to expect to be able to clear positions. Molecules allow positions to be "un-set" via removing conformers and Interchange.positions allows None to be passed. I'm an advocate against polymorphic inputs but I think this is a case in which it's intended. To keep the API cleaner, something line Topology.clear_positions could be introduced that does the same thing but ensures the array argument cannot be None. I'm happy to implement either solution.

Note also that get_positions allows None as an output:

In [1]: from openff.toolkit import Molecule

In [2]: assert Molecule.from_smiles("CCO").to_topology().get_positions() is None
# no error
j-wags commented

Thanks for the writeup. I'm in favor of allowing Topology.set_positions(None). Optionally Topology.clear_positions() would be fine in addition. Happy to take a PR for these :-)

mattwthompson commented

I decided for Topology.clear_positions() in #1827, which has the benefit of being less ambiguous about input types and the debatable downside of adding a new method. Since I'm the person who suggested Topology.set_positions(None) I'm also the person who can take it back