openff-toolkit

openff-toolkit copied to clipboard

Describe the bug Trying to create an openmm system with no nonbonded force results in an error due to the post-processing done here. In my case the nonbonded force is replaced by a custom nonbonded force implemented via smirnoff-plugins which creates the exclusions here before the toolkit so missing this step should be fine.

To Reproduce This is not quite the same as I have not added the custom force but the result is the same

from openff.toolkit.typing.engines.smirnoff import ForceField
from openff.toolkit.topology import Molecule

ff = ForceField("openff_unconstrained-1.3.0.offxml")
# remove all nonbonded handlers 
del ff._parameter_handlers["ToolkitAM1BCC"]
del ff._parameter_handlers["LibraryCharges"]
del ff._parameter_handlers["Electrostatics"]
del ff._parameter_handlers["vdW"]

water = Molecule.from_smiles("O")
system = ff.create_openmm_system(water.to_topology())

Output

IndexError Traceback (most recent call last) /var/folders/9q/nm__l0v13fggc72v94p0hybw0000gq/T/ipykernel_68062/4179821129.py in ----> 1 system = ff.create_openmm_system(water.to_topology())

~/miniconda3/envs/openff-plugins/lib/python3.8/site-packages/openff/toolkit/typing/engines/smirnoff/forcefield.py in create_openmm_system(self, topology, **kwargs) 1376 # TODO: Can we generalize this to allow for CustomNonbondedForce implementations too? 1377 forces = [system.getForce(i) for i in range(system.getNumForces())] -> 1378 nonbonded_force = [f for f in forces if type(f) == openmm.NonbondedForce][0] 1379 1380 nonbonded_force.createExceptionsFromBonds(

IndexError: list index out of range

Computing environment (please complete the following information):

• Operating system osx
• Output of running conda list
• openff-toolkit 0.10.0+32.gc316bfe0.dirty
• openff-toolkit-base 0.10.0
• smirnoff_plugins 0.0.1+4.g6118a88.dirty

Additional context

Thanks for opening this @jthorton. I'll keep thinking about a longer-term fix, but here's a quick monkey-patch solution using your reproducing example (just copying create_openmm_system and dropping the part where we expect a NonbondedForce to apply 1-4 scaling):

from openff.toolkit.typing.engines.smirnoff import ForceField
from openff.toolkit.topology import Molecule
import copy
from simtk import openmm

class MyForceField(ForceField):
    def create_openmm_system(self, topology, **kwargs):
        """Create an OpenMM System representing the interactions for the specified Topology with the current force field
        Parameters
        ----------
        topology : openff.toolkit.topology.Topology
            The ``Topology`` corresponding to the system to be parameterized
        charge_from_molecules : List[openff.toolkit.molecule.Molecule], optional. default =[]
            If specified, partial charges will be taken from the given molecules
            instead of being determined by the force field.
        partial_bond_orders_from_molecules : List[openff.toolkit.molecule.Molecule], optional. default=[]
            If specified, partial bond orders will be taken from the given molecules
            instead of being determined by the force field.
            **All** bonds on each molecule given must have ``fractional_bond_order`` specified.
            A `ValueError` will be raised if any bonds have ``fractional_bond_order=None``.
            Molecules in the topology not represented in this list will have fractional
            bond orders calculated using underlying toolkits as needed.
        return_topology : bool, optional. default=False
            If ``True``, return tuple of ``(system, topology)``, where
            ``topology`` is the processed topology. Default ``False``. This topology will have the
            final partial charges assigned on its reference_molecules attribute, as well as partial
            bond orders (if they were calculated).
        toolkit_registry : openff.toolkit.utils.toolkits.ToolkitRegistry, optional. default=GLOBAL_TOOLKIT_REGISTRY
            The toolkit registry to use for operations like conformer generation and
            partial charge assignment.
        Returns
        -------
        system : openmm.System
            The newly created OpenMM System corresponding to the specified ``topology``
        topology : openff.toolkit.topology.Topology, optional.
            If the `return_topology` keyword argument is used, this method will also return a Topology. This
            can be used to inspect the partial charges and partial bond orders assigned to the molecules
            during parameterization.
        """
        return_topology = kwargs.pop("return_topology", False)

        # Make a deep copy of the topology so we don't accidentally modify it
        topology = copy.deepcopy(topology)

        # set all fractional_bond_orders in topology to None
        for ref_mol in topology.reference_molecules:
            for bond in ref_mol.bonds:
                bond.fractional_bond_order = None

        # Set the topology aromaticity model to that used by the current force field
        # TODO: See openff-toolkit issue #206 for proposed implementation of aromaticity
        # topology.set_aromaticity_model(self._aromaticity_model)

        # Create an empty OpenMM System
        system = openmm.System()

        # Set periodic boundary conditions if specified
        if topology.box_vectors is not None:
            system.setDefaultPeriodicBoxVectors(*topology.box_vectors)

        # Add atom particles with appropriate masses
        # Virtual site particle creation is handled in the parameter handler
        # create_force call
        # This means that even though virtual sites may have been created via
        # the molecule API, an empty VirtualSites tag must exist in the FF
        for atom in topology.topology_atoms:
            system.addParticle(atom.atom.mass)

        # Determine the order in which to process ParameterHandler objects in order to satisfy dependencies
        parameter_handlers = self._resolve_parameter_handler_order()

        # Check if any kwargs have been provided that aren't handled by force Handlers
        # TODO: Delete this and kwargs from arguments above?
        known_kwargs = set()
        for parameter_handler in parameter_handlers:
            known_kwargs.update(parameter_handler.known_kwargs)
        unknown_kwargs = set(kwargs.keys()).difference(known_kwargs)
        if len(unknown_kwargs) > 0:
            msg = "The following keyword arguments to create_openmm_system() are not used by any registered force Handler: {}\n".format(
                unknown_kwargs
            )
            msg += "Known keyword arguments: {}".format(known_kwargs)
            raise ValueError(msg)

        # Add forces and parameters to the System
        for parameter_handler in parameter_handlers:
            parameter_handler.create_force(system, topology, **kwargs)

        # Let force Handlers do postprocessing
        for parameter_handler in parameter_handlers:
            parameter_handler.postprocess_system(system, topology, **kwargs)
        
        if return_topology:
            return (system, topology)
        else:
            return system


ff = MyForceField("openff_unconstrained-1.3.0.offxml")
# remove all nonbonded handlers 
del ff._parameter_handlers["ToolkitAM1BCC"]
del ff._parameter_handlers["LibraryCharges"]
del ff._parameter_handlers["Electrostatics"]
del ff._parameter_handlers["vdW"]

water = Molecule.from_smiles("O")
system = ff.create_openmm_system(water.to_topology())

The TODO two lines above seems quite relevant as a bit of debt 😅 .

Relates somewhat to discussions around revamping the SMIRNOFF spec, particularly https://github.com/openforcefield/standards/issues/2#issuecomment-930133340. Requiring at least one of NonbondedForce or CustomNonbondedForce seems to fit every use case I can think of. Any obvious pitfalls?

Requiring at least one of NonbondedForce or CustomNonbondedForce

I'm having trouble following here -- What object would have this requirement? The SMIRNOFF spec? Or the ForceField class/create_openmm_system method? I don't think it'd be a clean fit for either one.

Also, I don't think we'd need to require at least one -- Users should be free to have no nonbonded forces.

Generally, the section of code that this problem occurs in should not exist at all -- Each ParameterHandler should be able to unilaterally put its intended physics into the System. However the electrostatic and vdW terms are combined for performance reasons into the openmm.NonbondedForce. In a perfect world the vdWHandler and the ElectrostaticsHandler (+charge gen handlers) would deal separately with their own CustomNonbondedForces.

The simplest thing here may be to either:

• Refactor create_openmm_system into lower-level API points (maybe ForceField.create_forces, ForceField.postprocess_system, and ForceField.apply_exclusions) and expose those separately so that users like @jthorton can manually skip the step.
• Have a kwarg to create_openmm_system to skip applying the exclusions
• Automatically apply exclusions to any NonbondedForces or CustomNonbondedForces present in the system (this seems dangerous, since users like @jthorton may have already applied the ones they want)

I seem to be misremembering a detail I though was both in the spec and implementation. I agree that, though peculiar, force fields with no non-bonded interactions should be permissible.

I agree that, though peculiar, force fields with no non-bonded interactions should be permissible.

When computing gradients by finite difference I usually create a subset of a force field that only contains the relevant handler type. So if I was computing dU / d bond length then I would create a system from a force field only containing the BondHandler. Kind of a weird use case but definitely useful to be able to do.

I implemented this in Interchange; it was simpler than I expected and I only covered the simplest behavior in the tests, so there could be sharp edges when using it for interesting work.

This can be closed when the snippet above can run directly in the toolkit; I anticipate this will work after #1276 but it might take some small changes.

There are other issues with system subset creation but I think this one is resolved. I can run the original code snippet without error:

In [5]: from openff.toolkit.typing.engines.smirnoff import ForceField
   ...: from openff.toolkit.topology import Molecule
   ...:
   ...: ff = ForceField("openff_unconstrained-1.3.0.offxml")
   ...: # remove all nonbonded handlers
   ...: del ff._parameter_handlers["ToolkitAM1BCC"]
   ...: del ff._parameter_handlers["LibraryCharges"]
   ...: del ff._parameter_handlers["Electrostatics"]
   ...: del ff._parameter_handlers["vdW"]
   ...:
   ...: water = Molecule.from_smiles("O")
   ...: system = ff.create_openmm_system(water.to_topology())

/Users/mattthompson/software/openff-interchange/openff/interchange/smirnoff/_create.py:145: UserWarning: Automatically up-converting BondHandler from version 0.3 to 0.4. Consider manually upgrading this BondHandler (or <Bonds> section in an OFFXML file) to 0.4 or newer. For more details, see https://openforcefield.github.io/standards/standards/smirnoff/#bonds.
  _upconvert_bondhandler(force_field["Bonds"])